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SECTION  I 
INTRODUCTION 

In  recent  years,  the  development  of  mini  computers 
and  certain  algorithms  have  renewed  interest  that  a 
practical  way  may  be  found  to  extract  from  measured 
vibration  data  the  basic  structural  dynamic  properties 
which  govern  the  modal  response.  Such  an  ability  could 
greatly  enhance  the  usefulness  of  required  ground  vibration 
tests  of  new  or  modified  aircraft  for  evaluation  of 
aeroelastic,  aeroservoelastic , dynamic  loads  and  other 
dynamic  phenomena  directly  related  to  aircraft  safety. 

As  a consequence  of  these  developments,  a number  of 
methods  have  been  proposed  by  diverse  groups,  many  of 
which  are  outside  the  airframe  industry.  The  unfamiliarity 
of  some  of  the  methods,  semantics  problems,  and  proprietary 
considerations,  have  hindered  understanding  of  the 
basics  involved  for  evaluation,  development  or  adaptation 
for  more  specific  airframe  use.  This  report  describes  the 
foundations  of  some  of  the  more  prominent  methods  on  a 
common  basis  for  initial  comparisons  and  offers  a method 
based  on  the  steady-state  sinusoidal  response  rather  than 
the  transient  response  used  by  the  other  methods. 

We  assume  the  structure  is  modeled  by  a system  of 
second  order  differential  equations  with  constant  coefficients. 
In  matrix  and  vector  notation  this  system  of  equations  can 
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be  written  as 

Mx  ♦ Cx  ♦ Kx  ■ f 

The  matrices  M,  C and  K are  of  order  m and  are  called  the 

mass,  damping  and  stiffness  matrices  respectively.  The 

elements  of  these  matrices  are  real  numbers.  The  forcing 

function  or  excitation  f - f(t)  is  a vector  function  of  time 

with  m components  which  may  always  be  considered  as  known. 

The  response  x - x(t)  is  also  a vector  function  of  time  with 
i 

m components  which  are  determined  by  measurement.  The 
components  of  the  vector  functions  f(t)  and  x(t)  may  be 
complex  valued. 

] 

1 If  the  vector  function  u exp  (At)  satisfies  the  homogeneous 

. 

equation  Mx  ♦ Cx  ♦ Kx  * 0 then  the  complex  number  A is  called  a 
characteristic  value  and  the  vector  u a characteristic  or  modal 
vector  associated  with  A.  In  this  report  we  consider  methods 
for  determining  the  modal  parameters,  that  is,  the  characteristic 
values  A or  the  natural  frequencies  and  damping  coefficients 
and  modal  vectors  u.  The  natural  frequencies  (or  resonant 
frequencies)  and  damping  coefficients  (or  damping  ratios) 
are  readilv  obtained  from  the  characteristic  values  and 
conversely. 

The  methods  considered  in  this  report  have  the  capability 
of  determining  the  vibration  parameters  and  the  matrices  M, 

C and  K from  the  measured  responses  x(t)  to  excitations  f(t) 
which  have  only  one  (and  always  the  same)  component  different 
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from  zero.  For  most  of  these  methods,  only  symmetery  of  the 
matrices  M,  C and  K is  required.  It  is  not  necessary  to  assume 
any  relation  (e.g.  proportional  damping)  between  M,  C and  K. 

If  the  matrices  M,  C and  K are  not  symmetric  then, 
in  general,  one  must  excite  the  system  at  all  stations. 

That  is,  one  needs  a linearly  independent  set  of  excitation 
functions  f(t)  and  the  corresponding  responses  x(t). 

The  excitation  f(t)  and  the  resulting  response  x(t) 
are  experimental  quantities  determined  by  measurements. 
Accordingly  these  quantities  are  subject  to  error.  For 
purposes  of  this  report  we  assume  ideal  (error  free)  data. 

The  various  methods  were  examined  to  see  if,  at  least  in 
principle,  the  desired  quantities  could  be  obtained 
accurately  by  the  method.  The  sensitivity  of  a method  to 
errors  in  the  data  could  be  a decisive  factor  in  the  selection 
of  a method.  This  aspect  has  not  been  treated  because  the 
author  is  not  sufficiently  familiar  with  the  practical 
aspects  of  the  experimentation  procedures. 

In  Section  II  we  propose  a method  for  determining 
vibration  parameters  from  the  steady  state  response  to 
harmonic  excitations.  The  method  is  based  upon  F.q . (A36) 
of  Appendix  A.  We  show  how  to  compute  the  characteristic 
values  and  the  modal  vectors  u^,.  We  also  show  how  to 
cope  with  the  case  of  characteristic  values  of  multiplicity 
greater  than  1.  For  this  case  one  must  excite  the  system  at 
a number  of  different  stations  separately.  This  number  being 

the  same  as  the  multiplicity  of  the  characteristic  value. 
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The  method  of  Section  II  is  readily  extended  to  deal 
with  the  case  where  the  matrices  M,  C,  and  K are  not 
symmetric.  For  this  nonsymmetric  case  one  must  perform 
experiments  (separately)  at  each  station.  The  ability 
to  treat  the  nonsymmetric  case  seems  to  have  been 
overlooked  for  the  most  part. 

In  Section  III  we  describe  and  discuss  other  methods 
for  determining  vibration  parameters.  We  believe  that  we 
have  described  the  major  possible  methods  for  determining 
vibration  parameters  from  experiments  at  a single  degree 
of  freedom  when  no  relations  are  assumed  between  the 
matrices  M,  C,  and  K.  In  Section  IV  we  describe  and  give 
the  results  of  numerical  experiments  with  the  method  of 
Section  II.  We  note  possible  areas  for  further  investiga- 
tion and  present  our  conclusions. 

The  theoretical  treatment  may  be  carried  out  in 
either  the  time  domain  or  the  frequency  domain.  In 
general,  the  experimental  data  appears  in  the  time  domain. 
This  time  domain  data  is  transformed  to  frequency  domain 
data  by  numerically  Fourier  transforming  the  time  data. 
However,  if  one  measures  the  steady  state  response  to 
harmonic  excitations  then,  essentially,  one  obtains  the 
experimental  data  directly  in  the  frequency  domain.  Hence 
note  that  the  method  of  Section  II  can  also  be  used  in 
methods  which  determine  the  frequency  response  function 
experimental ly . 
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In  Appendix  A we  give  the  mathematical  background  for  the 
methods  described  in  this  report.  In  Appendix  B we  display  the 
connection  between  the  Fourier  Integral,  Fourier  Series  and  trig- 
onometric interpolating  polynomials.  Appendix  B is  background 
material  for  method  2 of  Section  III.  Appendix  C is 
background  material  for  methods  4A  and  4B  of  Section  III. 

Finally  in  Appendix  D we  give  instructions  for  the  user 
and  a program  for  the  method  of  Section  II. 


SECTION  II 


VIBRATION  PARAMETERS  FROM  THE  STEADY  STATE  RESPONSE 
TO  SINUSOIDAL  EXCITATIONS 

In  this  section  we  describe  a procedure  for  determining  the 
"complex  frequencies"  and  associated  complex  mode  shapes  of  a 
structure  from  its  measured  steady  state  response  to  sinusoidal 
excitations  of  frequencies  close  to  resonance.  In  principle, 
only  one  "point"  of  the  structure  need  be  excited  at  appropriate 
frequencies,  provided  there  are  no  multiple  characteristic  values. 
The  case  of  no  multiple  characteristic  values  is  treated  first 
after  which  we  show  how  to  deal  with  multiple  characteristic  values. 

Consider  the  following  system  of  linear  second  order 
differential  equations  wtih  constant  coefficients 

Mx  ♦ Cx  + Kx  = f (1) 

In  Fq.  (1)  M,  C and  K are  real,  square  matrices  of  order  m.  The 
matrices  M,  C and  K are  referred  to  as  the  mass,  damping  and 
stiffness  matrices  respectively.  The  vectors  x = x(t)  and  f = f(t) 
are  of  dimension  m with  time  dependent  components  which  may  be 
complex  valued. 

In  Appendix  A it  is  shown  that  the  steady  state  response  to  a 
harmonic  excitation  f(t)  = r exp  (iut)  of  a system  modeled  by  Fq. 

(1)  is  y exp  (iut)  where 

n j 

y = y(u,r)  = I u.v.r/(ia)  - A.)  (2) 

k»l  K K K 

For  k = 1 to  n = 2m  the  vector  functions  x^(t)  = u^  exp  (A^t)  are 
linearly  independent  and  satisfy  the  homogeneous  equation 

Mx  ♦ Cx  + Kx  = 0 (3) 
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Similarly,  the  vectors  z^Ct)  = exp  (A^t)  are  linearly  independent 
and  satisfy  the  transposed  homogeneous  equation.  The  vectors  u^ 
and  v^  satisfy  the  normalization  conditions  expressed  by  Eq.  (A17) 
of  Appendix  A.  The  complex  numbers  satisfy  the  characteristic 
equation 

det  [A2!^  + AC  + K]  = 0 (4) 


The  Eq.  (2)  is  valid  under  very  general  circumstances.  The  only 
requirement  for  the  validity  of  Eq.  (2)  is  that  for  a root  A^  of 

+ A^C  + K must  be  of  rank  m - m^. 

The  complex  numbers  A^  and  the  associated  vectors  u^  are  the 
"complex  frequencies"  and  mode  shapes  which  we  want  to  determine. 

The  procedure  to  be  described  below  extracts  these  quantities  from 
Eq.  (2).  For  the  nonsymmetric  case  the  vectors  v^  are  not  scalar 
multiples  of  the  corresponding  u^.  For  the  determination  of  the  v^. 
in  the  nonsymmetric  case  the  system  must  be  excited  at  each  degree 
of  freedom  with  frequencies  close  to  the  resonant  frequencies. 

In  principle,  the  procedure  described  below  is  general.  In 
the  application  of  interest  here  the  procedure  is  limited  to  systems 
for  which  the  characteristic  values  A^  lie  close  to  the  imaginary  axis. 
More  precisely  the  iteration  procedure  suggested  below  requires  values 


multiplicity  m^  the  matrix 


satisfying  the  relation 

xj  * V 


\ ' iwk  I < I ‘ ^ I for  a11 


For  the  case  in  which  the  matrices  M,  C and  K are  real  symmetric 


matrices, and  the  roots  A^  of  the  determinantal  equation,  Eq.  (4),  are 
single  and  well  separated,  the  vectors  v^  differ  from  the  u^  by  a 
numerical  factor,  at  most.  This  case  is  considered  now. 
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lot  ij  denote  the  in  dimensional  unit  vector  whose  i - component  is  1 
and  all  others  are  zero.  For  k * 1 to  in  and  ■ 1,2  let  v,  , exp  (im.  ,t) 
he  the  steads  state  response  to  i'j  exp  (iw^t).  .Accordingly  one  makes 
two  measurements  tor  each  characteristic  value.  The  vectors  are 
computed  from  experimental  data,  that  is,  from  the  response  to  the 
sinusoidal  excitation  i ^ sin  (m.  , t ) , see  Appendix  A,  Fq.  (\10). 

The  values  of  are  chosen  so  that 

u’kl  " Im  < ulkj  (S) 

and  , is  much  less  than  j.  (The  are  indexed  as  described 

in  Appendix  A.)  For  definiteness  we  discuss  the  case  where  the 
excitation  is  in  the  first  coordinate  direction  ij.  The  procedure 
is  the  samt'  regardless  of  the  fixed  coordinate  direction  excited. 

Here  now  we  describe  the  iterative  procedure  for  determining 
the  characteristic  values  and  the  associated  vectors  u^  and  v^. 

First  we  give  the  equations  for  computing  initial  values  for  . 
u^  and  Vj,.  Then  we  give  the  equat ions  for  computing  successive 
refinements  of  these  quantities  until  preassigned  tolerances  are 
attained. 

For  i * 1 to  m and  { = 1,  2,  we  have, on  neglecting  all  but 
one  term  from  Fq.  ( 2) , 

y j t " vlrl/(iu)j«  ' V (6) 


as  a first  approximation.  Set 


e - i (yj,  vn/vj, 


“u  ‘ sj 

T“.V  ' V 


where  the  sign  is  chosen  so  that  Im  |f.)  s 0.  (If  one  assumes  a value 
for  \j  and  computes  f.  from  Fq.  (7)  it  is  clear  that  this  is  the 
appropriate  choice  for  the  sign.)  From  Fq.  (7)  one  obtains 


(8) 


Aj  - i (£u>j  j -o)j2)/(C  - 1)  (8) 

for  j - 1 to  m.  Thus  Eq.  (8)  gives  a first  approximation  to  A ^ . This 
first  approximation  of  A . is  needed  for  computing  the  first  approximation 
of  the  vectors  u^  and  Vj.  From  Fxj.  (6)  we  obtain 

Yjl  - yj2  " ujVj’rli(“j2  ‘ ‘ XiUit°j2  ' V (9) 

Now  VjTj  » Vj.  is  the  first  component  of  the  vector  v^.  We  assume  the 
vector  Uj  normalized  so  that  its  first  component  v^.  = 1.  Then  from  F.q.  (9) 
we  have 

vIj  * <yljl  ‘ ylj2)(i“jl  - V<i»j2  - Vyi<“j2  • “jP  OH 
and  the  vector 

uj  ’ (yjl  * yj2)(iwjl  * V(iwj2  • Xj)/i(uj2  • wjllvl.i  (11) 

Hence  v.,as  a scalar  multiple  of  u.,is  given  by 

vj  ■ vuuj  (12) 

Thus  in  the  first  stage  the  computations  performed  are  those 
indicated  by  Fqs.  (7),  (81,  (10),  (11)  and  (12). 

Let  A1?,  u”,  and  vV  denote  the  new  value  which  one 
is  in  the  process  of  computing  and  let  A^,  Uj,  and  v^  denote  the  present 
value  of  these  quantities  respectively.  For  j = 1 to  m,  set 


yj«  ' yjt  Wl/(i“'.il  - V 


Compute 


'‘'P  A A A 1 / “> 

C ■ * Iyjlyil/yi2yi2l  * 


where  the  sign  is  chosen  as  in  Eq.  (7),  and 
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j 

L ■ 


I Xj 

I I 

I • 

I \ 


Then 


x1}  - i(^n  - wj2)/(£  - 1) 

■ 'iji'bi  ■ y 
u"  ■ yji«“n  ■ y'y 


(IS) 


(16) 

(17) 


and 


v? 

J 


v.  .u. 
1J  J 


(18) 


where  Aj  in  Eqs.  (16)  and  (17)  is  the  a”  computed  in  Eq.  (15), 
in  Eqs.  (17)  and  (18)  is  the  computed  in  Eq.  (16),  and  in  Eq.  (18)  u. 

is  the  u*|  computed  in  Eq.  (17).  At  the  same  time  one  is  computing 
A1?,  u"  and  vn  one  should  record  also 

J J J 


Cj 

Cj 


A1? 


(19) 

(20) 


and 


v"  x ^ 

Vj  j 


(21) 


These  identifications  should  be  made  in  the  first  stage  also. 

The  confutations  in  the  refinement  process,  Eqs.  (13)  through  (21), 
are  repeated  until  either  preassigned  tolerance  requirements  are 
satisfied  or  the  maximum  number  of  passes  through  the  refinement  process 
is  attained.  In  the  latter  event  one  is  faced  with  the  task  of 
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determining  why  the  tolerance  requirements  were  not  satisfied.  If 
the  tolerance  requirements  are  satisfied,  then  the  matrices  M,C,  and 
K are  determined  next. 

Before  discussing  the  computation  of  the  matrices  M,C,  and  K let 

us  consider  a modification  of  the  above  procedure  which  enables  us  to 

handle  the  case  of  multiple  characteristic  values,  for  definiteness 

suppose  x , ■ v j - x^  and  the  remaining  characteristic  values  are 

♦ 

simple,  let  v^  and  u. ^ denote  the  j—  component  of  the  vectors  Vj.  and 
Uj.  respectively.  According  to  the  discussion  in  Appendix  A we  may 
suppose 


V ■ V ■ v ■ v 

V21  12  v43  v34 


U11  " u12  " u33  " U44  “ 1 

As  for  the  case  of  no  repented  characteristic  values  the  remaining 
u^  vectors  are  assumed  normalized  so  that  the  first  component  has  the 
value  1. 

let  y(u,r)exp(ib>t)  denote  the  steady  state  response  to  the  harmonic 
excitation  rexp(iwt).  The  vector  y(w,r)  is  given  by  fiq.  (2).  Suppose 
we  have  obtained  experimentally  the  vectors  y(u>j.{,  r^)  for  * ■ 1»  2, 
see  Fq.  (15)  for  k ■ 1,  5,  6,  • • • , m and  also  the  vectors 
y(u'i4,r,),  y(u>jt,r^)  and  (y(io^{,r^).  Fran  this  data  one  obtains  a 
first  approximation  to  the  characteristic  values  x^  and  the  vectors  Uj.  for 
k - 1 to  m in  exactly  the  same  way  as  for  the  case  of  no  repeated  charac- 
teristic values. 

Thus  for  first  approx imat ions  compute  as  in  Hq.  (7) 

Cj  ■ (y  (<i,j j » (u,j j » )/y  (23) 


where  k ■ 1 for  j • 1,  5,  6,  • • • , m and  k - 3 for  j - 3. 

Then  a first  approximation  to  Xj  is  given  by 

Aj  " i(SWjl  ' a,j2)/(Cj  * 1}  (24] 

for  j » 1,  3,  5,  6,  •••,  m. 

Next  as  in  Eq.  (9)  we  have 

y("jl,rk)'y(“j2,rk)  ■ ujvj,rki(“j2  ' V(1“j2  ' V 

th  W- 

Let  >p (wje>rk)  denote  the  p — component  of  the  vector  y(Wj£,rk) • 

Then  from  Bq.  (25)  we  have 


vkj  ' Iyk("jl'rk)'ykl“j2’rk,;(1“jl  ' Vfl"j2  ‘ - ujj) 

(26) 

where  the  value  of  k mast  be  appropriately  assigned  depending  on 

j,  see  below  Bq.  (27).  Bq.  (26)  gives  the  k—  component  of  the  vector 

v. . 

J 

The  modal  vectors  are  obtained  as  before,  Eq.  (11).  Thus 

Uj  " [y (wj i » r^) *y (“)j 2 » rj() ] ( i-wj j " ■ Xj)/i(u'^2  ■“ii’vkj  <27' 

where  k • 1 for  j ■ 1,  5,  6,  •••,  m,  k ■ 2 for  j * 2 and  k = j for 
j - 3,  4. 

For  values  of  j for  which  X^  is  not  a multiple  characteristic 
value  the  vectors  Vj  are  given  by 

Vj  • Vj  (28) 

as  before.  For  values  of  j for  which  Xj  is  a multiple  characteristic 
value  the  associated  vectors  Vj  are  linear  combinations  of  the 
associated  Uj  vectors.  Rather  than  use  subscripted  subscripts  we 
chose  to  illustrate  the  procedure  for  determining  the  appropriate 
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b 


TOPI  11  > 


linear  combinations  by  example. 

For  the  example  being  described  the  vectors  Vj  and  v2  are  linear 


combinations  of  the  vectors  u and  u,. 

1 2 

Thus 


V1  “ all  “l  " a12  u2 1 
v2  " a21  U1  + a22  u2 1 


(29) 


Now  the  vectors  and  u2  are  known  from  Eq.  (27).  The  component  v^ 
of  Vj  and  the  component  v22  of  v2  is  known  from  Eq.  (26).  In  addition, 
because  of  the  special  fom  that  we  were  able  to  assume  for  the  vectors 
Vj  and  v2  we  have  that  the  com>onents  v21  ■ v12  - 0 for  the  vectors 
Vj  and  v2  respectively.  Utilizing  this  information  we  can  extract  a 
matrix  equation  from  the  system  of  vector  equations,  Eq.  (29)  for 
determining  the  coefficients  a^. . Thus  we  have 


u 


11 


u. 


12 


u 


21 


u. 


22 


‘11 


a 


21 


a 


12 


22 


'11 


0 


22 


(30) 


After  the  coefficients  a^  are  determined,  the  vectors  v^  and  v2  are 
completely  specified  (that  is,  the  first  approximation  thereof)  by 
Eq.  (29). 

In  a similar  fashion  one  has 


v5  ' "'ll  UJ  * b12  M 
v4  ' h21  u3  * b22  u4  I 


(31) 


From  Eqs.  (31)  one  obtains 


'33  u34 

bll  b21 

'43  u44 

b12  b22. 

33 


44J 


(32) 
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and  hence  the  vectors  and  are  completely  determined  also.  From 
the  above  discussion  it  should  be  clear  how  to  deal  with  any  number 
of  multiple  characteristic  values  whatever  their  multiplicity. 

The  computations  described  by  Eqs.  (23)  through  (32)  provide 
first  approximations  of  the  desired  quantities.  The  refinement  process 
for  the  case  of  multiple  characteristic  values  is  basically  the  same  as  for 
the  case  of  siirple  characteristic  values.  We  will  abide  by  the  same 
conventions  as  adopted  for  the  case  of  all  simple  characteristic  values 
see  Eqs.  (13)  through  (21). 

Set  as  in  Eq.  (13) 


ybrV  = y^e’V  - V 

A 

Then  replace  y(w-£,rk)  for  i - 1,  2 by  y(u>.£,rk)  in  Eq.  (23) 


(33) 


the  value  of  thus  obtained  compute 


With 


X"  ' ' V'Kj  * » 


(34) 


Then 


''kj " [y1^il.rk)Kicoil  - x.) 

uj = y^rrk)/vkj 


(35) 


(36) 


and  for  values  of  j for  which  X^  is  not  a repeated  characteristic  value 


“j  ‘ Vj 


(37) 


The  v^  corresponding  to  a multiple  characteristic  value  X ^ are  computed 
as  above,  see  Eqs.  (29),  (30),  (31)  and  (32). 

The  refinement  process  described  by  Eqs.  (33)  thru  (37) 
is  repeated  until  the  tolerance  requirements  are  satisfied. 
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Oice  the  characteristic  values  and  the  appropriately  normalized 
characteristic  vectors  are  known  the  system  matrices  M,C,  and  K are 
readily  obtained.  Set 

P(w)  = - w^M+iuC+K  (38) 

and 

Q(<*>)  * IV'fc/  (iw-Afc)  (39) 

Ne  noted,  Eq.  (A43)  of  Appendix  A,  that 

P(u.)  = Q'V)  (40) 

so  that 

P(W)Q(W)  ■ I (41) 

Let  the  prime  (')  denote  differentiation  with  respect  to  m. 

From  Eq.  (38)  we  have 

K = P(0) 

iC  » P'(0)  (42) 

-2M  - P"(0) 

From  Eq.  (41)  we  obtain 

P(0)  = Q_1 (0) 

P'(0)  -Q'^OiQ'WQ'^O) 

P”(0)  = -(2P'(0)Q'(0)+P(0)Q"(0))Q_1(0) 


(43) 


SECTION  III 

OTHER  METHODS  FOR  DETERMINING  THE  VIBRATION  PARAMETERS 

In  this  section  we  describe  and  discuss  other  methods  for  deter- 
mining the  vibration  parameters.  We  make  no  attempt  nor  claim  to 
discuss  all  methods.  Indeed  in  the  discussion  below  other  ways  to 
achieve  the  same  end  will  come  to  mind.  Each  variant  of  the  methods 
described  below  may,  of  course,  be  considered  another  method. 

The  methods  of  primary  interest  to  us  have  two  distinctive 
features;  First,  they  assume  no  relation  whatsoever  between  the  matrices 
M,  C,  and  K.  Secondly,  these  methdds  have  the  capability  of  determin- 
ing (almost)  all  vibration  parameters  of  interest  from  the  data 
arising  from  exciting  the  system  at  a single  degree  of  freedom,  at 
least  for  the  case  where  all  the  characteristics  values  are  different 
from  one  another. 

The  first  method  which  we  describe  requires  that  the 
system  be  excited  at  each  degree  of  freedom  individually.  Hence,  it 
is  not  a method  of  primary  interest.  However,  given  the  problem  of 
determining  the  matrices  M,  C,  and  K this  first  method  seems  to  be 
an  obvious  approach  to  the  solution.  It  shews,  at  least  in  principle, 
that  the  matrices  M,  C,  and  K can  be  determined  from  experimental  data. 

Let  us  excite  the  system  of  equations  Mx  + Cx  + Kx  = f at  the  kth 
coordinate  with  exp  (iw^t).  The  steady  state  response  will  be 

>k  exp  (iuijt)  for  k * 1 to  m.  (See  the  discussion  and  Eqs  (A35)  - (A40) 
of  Appendix  A . From  this  set  of  experiments  one  obtains  the  matrix 
equation 

-“i  M ♦ K + i-i  c = [yr..ymf  J (44) 
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The  Eq.  (44)  is  solvable  for  the  matrix  C.  We  have 


i 


C - (l/a^)Im{y1...yk)'J 


(45) 


If  this  set  of  experiments  is  repeated  with  an  i1  uij  then  we 
obtain  a system  of  matrix  equations 


m ♦ k - Re[yr..ym]'} 

’4  M+K  ’ Retyj—yJ^ 


(46) 


which  is  readily  solved  for  the  matrices  M and  K. 

Now  let  us  note  what  is  involved  in  this  first  method.  The 
experimental  requirement  is  a "shaker"  set  up  for  each  degree  of 
freedom  of  the  system  and  an  excitation  at  each  degree  of  freedom  of 
the  system  at  frequencies  uij  and 

The  computations  involved  in  this  method  are  the  determination  of 
the  colunns  of  the  matrices  [y^ - . -yml k > and  [y^. . -yml 2 (See  Eq*  (A40); 
computing  the  inverses  of  the  complex  matrices  (y^ . . . ym Jj»  and  . . -y  ] 
obtaining  the  matrices  M and  K (Eq.  46)  and  the  matrix  C (Eq.  (45)).  In 
addition,  if  one  needs  the  modal  vectors,  natural  frequencies  and  damp- 
ing ratios  then  one  has  to  perform  an  eigenvalue  and  eigenvector  analy- 
sis involving  the  matrices  M,  C,  and  K.  Stated  in  another  way  one  has 
to  determine  a linearly  independent  set  of  n functions  of  the  form 
u exp  (At)  which  satisfy  the  equation  Mx  ♦ Cx  ♦ Kx  ■ 0. 

In  the  determination  of  the  matrices  M,  C,  and  K there  are  no 
simplifying  assumptions.  None  of  the  operations  or  computations  are 
performed  approximately.  Indeed  even  symmetry  of  the  matrices  M,C, 
ami  K is  not  needed  or  used. 
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The  value  of  the  components  of  the  complex  vectors  depend  upon 
the  frequency  u>  . As  usual,  let  y^  denote  the  jth  component  of  the  vector 
y^.  if  one  had  a graph  of  jy.jJ  as  a function  of  u>  one  would  observe 
that  |y.jj  is  small  when  u>  is  not  close  to  a natural  frequency  of  the 
system.  This  remark  is  evident  from  Fq  (A36) . Hence  the  inverses 
[y  ...y  ] j and  [ y ^ . . . ym ] *1  obtained  from  experimental  data  may  be  very 
inaccurate  for  most  choices  of  u>  and  u> 

In  the  next  method  considered  the  frequency  response  function  is  de- 
termined experimentally  as  a function  of  frequency  u>,  [l-b].  This  method  is 
characterized  by  its  use  of  techniques  of  Fourier  Analysis.  An  analo- 
gous method  based  on  the  Laplace  transform  is  also  possible. 

LettF  f(t))  denote  the  Fourier  transform  of  the  function  f(t).  We 
have  [u>  2 M ♦ iu>C  + K]F<x(t))  = F(f(t)>or 

Fix(t)}  = [-u>2  M + iwC  + K]"1  Flf  (t)  > (47: 

Now  let  fj(t)  denote  a vector  valued  function  of  t with  all  but  the 

jth  component  identically  zero.  Also,  let  x^  (t)  denote  the  response 

2 -1 

to  fj(t).  Then  the  jth  colunn  of  [-w  M + itoC  + K]  is  given  by 


2 -1 
jth  column  of  [-w  M ♦ i u>C  ♦ K] 


F{\j  (t) } / F{ fj (t) ) (48) 


2 -1 

For  j running  through  the  values  from  1 to  m the  matrix  [-u“  M + iwC  + K] 

is  completely  specified  as  a function  of  u , see  Eq  (A43). 

If  the  matrices  M,  C,  and  K are  synmetric  and  if  the  natural 

frequencies  of  the  system  are  well  separated  then  the  vibration  parameters 

2 -1 

can  be  determined  from  a single  column  of  [-w  M + iwC  + K]  such  as  given 


i 


by  Eq  (48).  For  example,  from  Eq  (A43)  the  first  colunn  of 
2 - 1 

[-W  M + i (>♦  K]  is  the  expression 


(49) 


Thus  the  vectors  u^v^  and  the  complex  numbers  ^ for  k = 1 to  n 
need  to  be  determined  so  that  the  equality 


l Viv/Ci-V  = F{xi  Ct)  >/F(  (t) ) (50) 

k=l 

holds  (at  least  in  a curve  fitting  sense). 

This  second  method  is  suited  to  an  impulse  type  excitation  but 
other  forms  of  excitation  may  be  used  (4  ].  The  experimental  pro- 
cedures for  these  first  two  methods  differ.  In  this  second  method  the 
transient  part  of  the  response  is  inportant  in  the  analysis  process, 
whereas,  in  the  first  method  the  steady  state  part  was  used.  Note 
that  the  method  of  Section  II  uses  the  same  kind  of  experimental  data 
as  the  first  method  of  this  section. 

For  this  second  method  we  observe  that  there  are  two  computational 
tasks.  First,  a colunn  of  the  frequency  response  function  must  be 
computed,  Eq  (48),  from  the  excitation  and  response  data.  Secondly, 
the  modal  vectors  and  characteristic  values  must  be  extracted  from  the 
computed  frequency  response  function  data.  Let  us  examine  these  two 
tasks  further. 

For  the  first  task  the  excitation  data  and  the  components  of  the 
vector  response  data  are  fitted  by  trigonometric  interpolation  poly- 
nomials to  obtain  the  frequency  response  function  data,  Appendix  B. 

If  the  system  of  differential  equations  Mx  ♦ Cx  ♦ Kx  * f is  of  order 
m then  one  has  up  to  m+1  fits  by  trigonometric  polynomials,  depending 
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upon  the  nunher  of  degrees  of  freedom  of  interest.  l;or  a complete 
analysis  (determination  of  the  matrices  M,  C,  ;uid  k)  all  m+1  fits  must 
he  made.  It  is  recognized  that  the  frequency  response  function  data 
is  the  result  of  certain  approximations. 

A number  of  ways  for  extracting  the  characteristic  values  and 
vectors  from  the  frequency  response  function  data  are  noted  in  |j*4] . 
l'he  problem  of  course  is  the  determination  of  the  quantities  u^v^  and 

so  that  the  expression  (49)  fits  the  frequency  response  function 
data.  Present  practice  is  to  least  squares  fit  the  frequency  response 
function  data  by  the  expression  (49).  The  "normal"  equations,  result- 
ing from  the  direct  approach  to  a least  squares  fit  of  the  expression 
(49)  to  the  frequency  response  function  data,  are  not  solvable  exactly. 
Thus  the  modal  vectors,  natural  frequencies  anti  damping  ratios  are 
determined  approximately.  The  graph  of  the  expression  (49),  after  the 
quantities  UjA  lh  :UK*  *iave  ^een  determined,  when  compared  with  the 
graph  of  the  frequency  response  function  data,  gives  a visual  check  of 
the  fitting  process. 

It  is  clear  from  the  above  discussion  that  there  is  considerable 
data  manipulation  and  processing  associated  with  this  second  method. 

'[he  desired  vibration  parameters  are  obtained  by  numerical  (approximate) 
methods.  On  the  other  hand  this  second  method  may  not  be  as  sensitive 
to  errors  in  the  data  as  some  of  the  other  methods. 

In  method  2 (just  considered)  the  vibration  parameters  are 
obtained  from  the  measured  (experimentally  determined)  frequency  response 
function.  In  method  3 to  be  considered  next,  the  vibration  parameters 
are  obtained  from  the  measured  impulse  response  function.  This  method 
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is  based  on  Eq  (A31)  when  expressed  as  Eq  (A45),  Appendix  A. 

When  the  forcing  function  f(t)  has  only  one  component,  for 
example,  the  first  component  different  from  zero  then  Eq  (A45)  sim- 
plifies to 

t 

x(t)  = JQ  I1(t-T)f1(T)di  (51) 

Set  t^  - tjc  ■ kAT.  Observe  that  Ijftj-T^)  = I-^Ct ^ for  k^j.  If  the 
integral  in  Eq.  (51)  is  evaluated  by  the  trapezoidal  rule  we  have 

X(t.j)  = (Ijftj)  ^(0/2  + l I^x^)  ^(t^At  (52) 

k-1 

Note  that  for  j=l  Eq  (52)  is  solvable  for  1^ (t^) . Now  that  I^fx^) 
is  known,  then  for  j=2  Eq  (52)  is  solvable  for  I^^)  so  on*  In 
this  way  the  impulse  response  function  1-^  (t)  is  determined  from  experi- 
mental data. 

The  impulse  response  function  data  is  used  to  determine  the  vibra- 
tion parameters.  From  Appendix  A,  Eq  (A34)  we  have 

n 

Il(t)  = E ukvikexp^kt)  (53) 

k 1 

The  quantities  u,  vlk  and  A^  need  to  be  determined  so  that  the  function 
defined  by  Eq  (53)  fits  the  impulse  response  function  data  obtained 
from  Eq(52). 

We  will  now  describe  briefly  a way  to  obtain  the  vibration  para- 
meters from  the  impulse  response  function  data.  This  procedure  is  given 
also  in  [7,  pp  270-280].  A component  of  Ij(t)  is  of  the  same  form  as 
I^(t)  as  given  by  Eq  (53).  That  is,  the  jth  component  of  the  vector 
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valued  function  Ij(t)  is  scalar  function  of  the  form 


S(t)  = l a.  (exp(A.t) 


k=l 


(54) 


Briefly  stated,  the  problem  is  to  determine  the  complex  numbers  Aj. 
and  the  coefficients  a^  from  the  data  t . = t+jh  and  g.  = g(tj)  for 
j = 0,1,2  and  so  on. 

The  problem  is  solved  in  two  stages.  In  the  first  stage  the 
Aj.'s  are  determined  and  in  the  second  stage  the  coefficients  a^. 
First,  we  determine  the  coefficients  C,  of  a difference  equation 

g(t)  ♦ Cjglt+h)  >...+  Cng(t+nh)  = 0 (55) 


which  the  functions  exp  (A^t)  (and  hence  any  linear  combination  of  the 
exp  (A^t)  also)  will  satisfy  for  a fixed  value  of  h and  all  values  of 
t.  Using  the  data  we  obtain  the  system  of  equations 

llRk  * C2Rk+l  lnpk+n-l  = pk-l  (5b) 

for  k=l,  ...  , n which  is  solved  for  the  C.  s. 

.1 

Now  set 


then 


p = exp  (Ah) 


exp  (Ajh)  = p-' 

IVe  see  that  g(t)  = exp  (At)  will  satisfy  Fq.  (55)  if  p is  a root  of 
the  polynomial  equation 

Cnpn  +...♦  ciP+1  * 0 (57) 
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Thou 


\h  - logep 


(SK) 


Thus  the  \y  are  determined  by  the  n roots  p^  of  Eq.  (57). 

The  n roots  pj.  of  Fq.  (57)  ;ire  used  also  in  the  determination 
of  the  coefficients  a^.  From  Fq.  (54)  we  obtain  the  system  of 
equations 


Vi 


k 

a p 
ir  n 


(Sh) 


for  k ■ 0,1,  ....  n-l.  Iliis  system  of  equations  determines  the 
coefficient  a^.  A method  for  computing  the  values  of  the  a^’s  from 
Hqs.  (59)  is  given  in  [7,  p 274 1. 

Let  p^  denote  an  n dimensional  vector  whose  components  are  l, 

p^ p”  * respectively,  where  p^  is  a root  of  Fq.  (57).  Ihe 

coefficient  matrix  for  the  l^s,  Fq.  (5(0  can  he  written  in  the 
dyadic  form 


J,  Vk CX|’  (Y>  i\  i’i 


((>0) 


It  is  clear  from  this  expression  that  the  coefficient  matrix  for  the 
t^s  is  nonsingular  provided  the  roots  p^  of  Eq.  (57)  are  all 
different.  This  expression  may  also  he  helpful  in  examining  the 
coefficient  matrix  of  the  C^s  for  ill  conditioning. 

The  fact  that  the  characteristic  values  \,  are  all  different  does 
not  guarantee  that  the  coefficient  matrix  for  the  f^s  is  nonsingular. 
For  example,  if  the  characteristic  values  A.  and  A^  and  the  sample 
interval  h satisfy  the  condition 
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((>1) 


h(Ak-A .)  = 2n i 

then 

Pj  = oxp(X.h)  = exp(XjJi)  = pk  (62) 

.iikI  the  coefficient  matrix,  by  hq.  (no),  is  singular. 

let  review  briefly  the  computations  involved  in  this  method, 
first,  there  are  the  numerical  integrations  in  the  deterininat  ion  of  the 
impulse  response  function  data.  Secondly,  the  linear  vetom  riven  hv 
l.q.  (r-h)  must  be  solved  loi  the  difference  equation  coeflicients 
| rrors  in  the  coefficient  matrix  and  in  the  rirlit  hand  side  of 
fq.  (5t» ) will  affect  the  l.-'s.  Thirdly,  the  roots  of  I <t . (57) 
must  be  determined.  Lastly,  the  linear  system  pi  von  hv  lq.  ( S*-» ) 
must  be  solved  for  the  amplitudes  for  each  decree  ol  freedom 
of  interest. 

We  will  refer  to  the  next  two  methods  as  method  L\  and  III 
respectively,  [S,  P,  1 0 1 . In  these  methods  the  initial  (or 
original)  data  is  a solution  x(t)  of  the  homogeneous  equ.it  ion 
Mx  » Cx  + Kx  = 0.  I he  vector  function  x(t)  is  a linear  combination 
of  the  functions  u^explX^t).  Huts  x(t)  is  ot  the  same  form  as 
the  impulse  response  function. 

It  is  clear  that  the  procedures  described  in  method  3 could  lx' 
used  to  determine  the  characteristic  values  and  associated  charac- 
teristic vectors  from  the  function  x(t).  The  function  x(t)  is  original 
experimental  data  whereas  the  impulse  response  function  is  computed  from 
experimental  data.  Thus  the  numerical  integration  step  of  method  3 used 
to  obtain  the  impulse  response  function  could  be  eliminated  in  the  case 
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where  x(t)  is  a free  response  of  the  structure.  However,  the  system 
matrices  M,  C,  and  K are  not  computable  from  the  results  so  obtained  (see 
remarks  which  follow  Eq  (A43)  of  Appendix  A) . 

The  equations  and  details  of  computations  in  the  methods  4A 
and  4B  are  given  in  Appendix  C.  Here  we  only  need  to  note  what  is 
involved  in  these  two  methods  and  contrast  and  compare  them  with  the 
other  methods  whenever  possible. 

The  first  step  in  Methods  4A  and  4B  is  the  determination  of  the 
"system”  equations.  The  second  step  is  the  determination  of  the 
characteristic  values  and  vectors  from  the  system  equations  obtained 
in  the  first  step.  All  the  previous  methods,  except  for  the  very 
first  of  this  section,  first  determined  the  characteristic  values  and 
a particular  normal ization  of  the  characteristic  vectors  so  that  the 
system  equations  can  be  determined  from  these  quantities. 

In  mctliod  4A  the  system  equations  determined  are  the  equations 
of  an  equivalent  first  order  system,  Fqs  (C4)  and  (C5)  of  Appendix  C. 
Thus,  the  matrices  M,  C,  and  K are  not  obtained  in  this  method.  A 
numoer  of  numerical  integrations,  see  Fqs. (C20)  and  (C21),  arc  needed 
to  produce  the  data  required  in  the  determination  of  the  system 
equations,  see  Fq  (C9) . 

The  method  4B  is  free  of  the  numerical  integrations  required 
in  method  4A.  This  freedom  is  achieved  by  determining  a system  of 
difference  equations  which  has  the  same  solutions  as  the  differen- 
tial equations  (see  the  discussion  following  Fq  (C22)  of  Appendix  C). 

The  matrices  M,  C,  and  K are  not  obtained  in  method  4B. 

After  the  system  equations  are  obtained  an  eigenvalue  and 


eigenvector  analysis  are  performed  to  obtain  the  characteristic 
values  and  vectors.  One  should  note  that  the  matrices  involved  in 
the  computations  in  methods  4A  and  4B  are  twice  as  large  as  those 
in  the  other  methods. 

In  the  preceding  pages  of  this  section  we  have  described  methods 
2,  3,  4A  and  4B.  These  methods  essentially  satisfy  our  self  imposed 
requirements.  That  is,  the  matrices  M,  C,  and  K are  not  assumed 
linearly  connected  and  the  vibration  parameters  can  be  determined  from 
the  response  to  the  excitation  at  a single  point.  In  the  discussion  of 
these  methods  one  acquires  some  awareness  of  computational  require- 
ments. However,  the  previous  discussion  was  primarily  concerned  with 
the  foundations  of  the  methods.  Here  now  we  want  to  make  quantitative 
estimates  of  some  computational  costs  associated  with  two  of  the  methods. 

The  Fast  Fourier  Transform  Algorithm  provides  a significant  reduc- 
tion in  the  cost  of  the  computations  indicated  by  Fq.  (48).  Thus  an 
operations  count  for  the  process  which  gives  the  value  of  a component 

of  Ffx- (t) } at  N stations  is  of  the  order  of  N log-N.  In  method  3 
J ^ 

the  determination  of  the  value  of  a component  of  I^(t)  at  N stations  is 
the  same  as  solving  a triangular  system  of  linear  equations.  Hence, 
the  operations  count  is  of  the  order  of  N(N+l)/2. 

The  method  of  Section  II  and  the  methods  2 and  3 can  determine 
the  natural  frequencies,  damping  ratios  and  corresponding  mode  shapes 
at  specified  degrees  of  freedom  of  interest.  That  is,  (we  haven't 
emphasized  this  point)  one  doesn't  have  to  compute  these  vibration 
parameters  for  all  m degrees  of  freedom.  However,  the  system  matrices 
M,  C,  and  K can  be  determined  by  the  method  given  above  only  if  the 
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components  of  the  appropriately  normalized  characteristic  vectors 
are  known  for  all  m degrees  of  freedom. 

In  methods  4A  and  4B  "equivalent"  system  matrices  of  order  n = 2m 
are  determined  first.  Thus,  one  must  select  m degrees  of  freedom  and 
collect  and  use  the  response  data  at  all  of  these  degrees  of  freedom. 

The  selection  of  the  data  collection  stations  and  the  sampling  frequency 
may  influence  the  calculation  of  the  characteristic  values  and  associated 
characteristic  vectors.  There  are  similar  difficulties  associated  with 
the  other  methods. 
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SECTION  IV 

RESULTS  AND  CONCLUSIONS 

In  this  section  we  describe  the  numerical  experiments  used  in 
testing  the  method  of  Section  II  and  give  the  results  of  these  experi- 
ments. We  note  again  some  of  the  observations  made  in  the  previous 
sections. 

Sample  problems  to  test  the  method  of  Section  II  are  readily  con- 
structed. One  uses  Eq.  (A36)  of  Appendix  A for  this  purpose.  Thus, 
in  Eq.  (A36)  we  assigned  appropriate  values  to  x^  and  u^  for  k=l  to 
m.  This  in  turn  fixes  the  values  of  x^,  uk  and  v^  for  k=l  to  n=2m. 
Then  for  a given  vector  r and  values  of  w Eq  (A36)  gives  vectors  y. 

We  used  the  vectors  y,  obtained  in  this  manner,  in  the  method  of 
Section  II  to  compute  the  x^,  uk  and  v^. 

For  a system  with  m=4  the  computed  values  of  x^,  and  v^  agreed 
with  the  assigned  values  to  about  8 decimal  places  after  10  or  so 
iterations.  In  Appendix  D we  display  the  computed  values  and  assigned 
values  in  one  case  where  the  characteristic  values  are  all  different 
and  in  a second  case  for  a characteristic  value  of  multiplicity  2. 
These  results  show  that  the  method  of  Section  II  is  computationally 
feasible  and  capable  of  high  accuracy. 

The  methods  discussed  in  Section  III  are  also  capable  of  high 
accuracy,  at  least  in  principle.  However,  as  is  well  known,  the 
excessive  cost  of  high  accuracy  frequently  results  in  compromise. 
Limitations  on  computational  capabilities  also  limit  the  attainable 
accuracy. 

The  method  of  Section  II  requires  the  steady  state  response 
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to  a sinusoidal  excitation.  Method  2 computes  the  frequency  response 
function  from  Impulse  or  random  excitation.  Method  3 computes  an 
Impulse  response  function  also  from  an  impulse-like  excitation.  Methods 
4A  and  4B  are  based  upon  the  free  response.  Thus,  the  way  in  which 
a structure  is  excited  influences  the  choice  of  method  for  computing 
vibration  parameters.  Alternately,  first  choosing  a method  for  com- 
puting vibrations  parameters  determines  how  the  structure  is  excited. 

From  these  remarks  it  is  clear  that  the  various  methods  comple- 
ment one  another.  We  have  noted  that  the  method  of  Section  II  assumes 
the  structure  is  lightly  damped.  If  this  is  not  the  case  then 
methods  3 or  4A  or  4B  are  possibly  better  suited  to  determining  the 
vibration  parameters.  On  the  other  hand  if  one  wants  the  matrices 
M,  C,  and  K then  a method  other  than  method  4A  or  4B  should  be  used. 

In  all  the  methods  there  is  some  stage  which  involves  the 
determination  of  the  characteristic  values  A^  and  the  modal  vectors 
u^.  In  the  method  of  Section  II,  the  numbers  A^  and  the  vectors 
^ are  determined  directly  from  the  experimental  data.  In  the 
other  methods  certain  intermediate  steps  are  necessary  before 
calculating  the  A^  and  u^.  Thus  in  the  methods  4A  and  4B  one  first 
determines  the  matrices  of  a system  equivalent  (in  the  sense  noted 
earlier)  to  the  system  given  by  the  matrices  M,  C and  K.  The  A^ 
and  u^  are  then  obtained  as  the  solution  to  the  generalized  algebraic 
eigenvalue  problem  defined  by  the  matrices  of  the  "equivalent"  system. 

In  method  3 one  first  determines  an  impulse  response  function 
which  we  denoted  by  Ij(t).  The  function  r (t)  is  defined  as  the 
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solution  of  a Volterra  integral  equation  of  the  first  kind,  Fq.  (51). 

By  transforming  to  the  frequency  domain,  the  integral  equation  problem 
of  method  3 is  reduced  to  the  algebraic  problem  of  method  2.  See 
[11].  The  impulse  response  function  (or  its  transform,  that  is,  the 
frequency  response  function, is  assumed  to  have  a specific  form, 

I!q.  (53)  or  Fq.  (50)).  After  either  the  impulse  response  or  the 
frequency  response  function  is  known  the  identification  of  parameters 
X^  and  u^  begins. 

Our  purpose  in  this  report  was  to  present  (Section  II)  a method 
for  determining  the  characteristic  values,  modal  vectors  and  the  mass, 
damping  and  stiffness  matrices  M,  C,  and  K from  the  steady  state  response  to 
sinusoidal  excitations  at  frequencies  close  to  resonance.  The  mathe- 
matical background  for  this  method  (Appendix  A)  with  a few  additions 
is  also  the  background  for  other  methods  which  we  described  and  discussed 
in  Section  III.  Thus,  we  were  able  to  note  difficulties  (e.g.,  an  inte- 
gration being  replaced  by  a numerical  process)  with  each  method. 

The  effects  of  these  difficulties  need  to  be  investigated  further. 

These  effects  could  be  determined  by  analysis  but,  possibly,  more  readily 
by  numerical  experiments.  Thus  nunerical  experiment  could  determine  the 
computer  resources  required  by  each  method  to  obtain  the  solution  to  a 
prescribed  accuracy.  Secondly,  the  sensitivity  of  each  method  to  error 
in  original  data  could  be  estimated  also  by  numerical  experiment. 


30 


APPENDIX  A 

SYSTEMS  OF  LINEAR  DIFFERENTIAL  EQUATIONS  WITH 
CONSTANT  COEFFICIENTS 

In  this  section  we  give  the  mathematical  background  (notation 
and  assumptions)  for  procedures  for  determining  mass,  stiffness  and 
dancing  properties  of  a structure.  The  mathematics  required  is  a 
knowledge  of  the  nature  and  properties  of  the  solutions  to  a system 
of  second  order  linear  differential  equations  (Lagrange's  equation) 
with  constant  real  coefficients.  The  determining  of  the  structural 
parameters,  mathematically  speaking,  is  of  the  nature  of  an  inverse 
problem.  Specifically,  we  are  concerned  with  determining  the 
differential  operator  (those  constant  coefficients).  Such  a system  of 
differential  equations  can  be  characterized  by  the  solutions  of  the 
homogeneous  equation  and  the  transposed  homogeneous  equation.  The  roots 
of  its  characteristic  equation  and  the  associated  characteristic  vectors 
are  required  for  one  such  characterization. 

Below  we  develop  expressions  for  solutions  of  a system  of 
second  order  linear  differential  equations  with  constant,  real 
coefficients.  These  expressions  are  obtained  from  expressions  for 
the  solutions  to  an  equivalent  first  order  system  of  equations.  The 
results  obtained  are  found  in  texts  which  discuss  systems  of  ordinary 
differential  equations  in  some  detail,  however  not  always  in  the 
explicit  form  given  below.  For  additional  details  concerning  systems 
of  differential  equations  see  references  [12,  13,  14 J 


[Mu  ♦ CD  ♦ K]  x (t)  - f(t) 
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where  D indicates  differentiation  with  respect  to  t,  denote  a system 
of  second  order  differential  equations.  The  mass,  damping  and 
stiffness  matrices,  M,  C and  K,  are  square  matrices  of  order  m with 
real  entries.  The  vectors  x * x(t)  and  f * f(t)  are  of  dimension 
m with  time  dependent  components  which  may  be  complex  valued.  The 
matrices  M,  C and  K are  usually  synmetric  but  this  condition  is  not 
essential  for  the  development  given  here. 

It  is  assuned  that  the  homogeneous  equation 

Mx  + Ck  + Kx  = 0 (A2) 

has  n = 2m  linearly  independent  solutions.  The  associated 
characteristic  equation  is 

det  [MA2  ♦ CA  ♦ K]  - 0 (A3) 

The  left  hand  side  of  Eq.  (A3)  is  a polynomial  in  A with  real 
coefficients  . This  polynomial  is  of  actual  degree  n.  Since  the 
coefficients  of  the  characteristic  equation  are  real,  the  complex 
conjugate  of  each  complex  root  of  Eq.  (A3)  is  also  a root. 

If  A is  a simple  root  of  the  characteristic  equation  then  there 
is  a solution  of  the  homogeneous  equation  of  the  form 

x(t)  = u exp  (At)  (A4) 

where  u is  a nontrivial  constant  vector.  The  vector  u satisfies  the 


condition 


[MA"  + CA  + K]  u = 0 


and  is  unique  except  for  a non- zero  multiplicative  factor.  If  A is 

2 

a root  of  multiplicity  k > 1 we  assume  the  matrix  MA  + CA  + K is  of 
rank  m - k.  In  this  event  there  is  a linear  manifold  of  vectors  of 
dimension  k and  every  vector  u in  this  linear  manifold  satisfies  Eq.  (A5) . 
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Accordingly  for  every  such  vector  u,  x(t)  ■ u exp  (Xt)  is  a solution 
of  the  homogeneous  equation. 

The  roots  of  the  characteristic  equation  are  oftentimes  referred 
to  as  characteristic  values  or  eigenvalues.  The  vectors  u which 
satisfy  Eq.  (A5) , will  be  referred  to  as  the  modal  or  characteristic 
vectors  or  eigenvectors  corresponding  to  the  characteristic  value  X. 

The  characteristic  equation 

det  [MTX2  + CTX  ♦ KT]  = 0 (A6) 

of  the  transposed  homogeneous  equation 

MTx  + CT*  + KTx  = 0 (A7) 

has  exactly  the  same  roots  as  Eq.  (A3).  For  a given  X the  matrices 
MX2  ♦ CX  ♦ K and  M*X2  ♦ CTX  + have  the  same  rank.  Consequently, 
the  dimensions  of  the  solution  sets  of  Eqs.  (A2)  and  (A7) , corresponding 
to  a characteristic  value  X,  are  the  same.  If  the  matrices  M,  C and  K 
are  symmetric  then  the  solution  sets  are  identical.  In  either 
event  there  are  n characteristic  values  Xk  and  associated  constant  vectors 
u^  and  vk  such  that  for  k = 1 to  n 

xk(f)  = \ exp  (Akt) 
and 


zk(t)  * vk  exp  (Xkt) 

are  linearly  independent  solutions  of  Eqs.  (A2)  and  (A7)  respectively. 
Set  Xj(t)  = x(t)  and  x2(t)  = *^(t).  Then  the  system  of  first 


order  equations 


■ 

I 

0 

X1 

♦ 
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-I 

X1 

0 
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• 

- - 

(A8) 
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is  equivalent  to  the  system  given  by  Eq.  (Al) . If  the  matrices  M,  C, 
and  K are  synmetrical  one  can  obtain  an  equivalent  first  order  system 
which  is  also  synmetrical  [15] . Express  Eq.  (A8)  symbolically  as 

A?  + By  = f (A9) 


Basic  or  fundamental  solutions  of  the  homogeneous  equation 
Ay  + By  = 0 


(A10) 


are  of  the  form  u exp  (At) , where  u is  a vector  of  dimension  n.  It 
is  clear  that  a linearly  independent  set  of  modal  vectors  is  given  by 


LVkJ  (A11) 

for  k = 1 to  n,  where  A^  and  u^  are  the  characteristic  values  and 
associated  modal  vectors  of  a set  of  linearly  independent  solutions  of  Eq 
(A2) . The  characteristic  values  Afc  satisfy  the  characteristic  equation 


det  [AA  + B]  = 0 

and  the  associated  modal  vectors  satisfy  the  condition 
[XkA  + B]u^  = 0 

For  the  transposed  homogeneous  first  order  system 

T T 
A V + B y = 0 

the  modal  vectors  v^>  f°r  k = 1 to  n,  are  given  by 


(A12) 


(A14) 


(Al  5) 


In  Eq.  (A15)  A^  and  v^  denote  the  characteristic  values  and  modal 
vectors  of  a set  of  linearly  independent  solutions  of  Eq.  (A?) . 

The  stiffness  matrix  K is  not  required  to  obtain  the  v^  from  v^. 

A modal  vector  v^  corresponding  to  the  characteristic  value  A^ 
satisfies  the  condition 

Xk  *kA  + vJb  = 0 (A16) 
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From  Eqs.  (A13)  and  (Alb)  it  follows  that  if  j + k and  X.  f X,, 

J * 

then  the  associated  modal  vectors  Vj  and  Uj  satisfy  the  orthogonality 
conditions 


J=AQk 


T 

0 and  VjBu^ 


(A17) 


F.ven  though  X^  takes  on  the  same  value  for  several  consecutive  values 
of  the  index  k,  the  associated  model  vectors  u^  and  v^  can  be 
determined  so  that  the  orthogonality  conditions  Eq.  (A17)  still  hold 
for  j i k.  The  vectors  u^  and  v^  may  be  normalized  so  that  in 
addition  to  the  orthogonality  conditions  they  satisfy  also  the  conditions 


vJa^  = 1 and  vJbl^  = -Xk  (A18) 

Write 

U = [Uj  •••  ur]  and  V = [v^  •••  v^]  (A19) 

Here  U denotes  a matrix  whose  kth  column  is  the  vector  u^. 

Similarly,  the  kth  column  of  the  matrix  V is  the  vector  v^.  Let 
A denote  a diagonal  matrix  with  the  characteristic  values  X^  in  the 
diagonal  positions.  The  relations  given  in  Eqs.  (A17)  and  (A18)  can 
be  expressed  in  matrix  form  also,  namely 


VTAU  = I and  VTBU  = -A 


(A20) 


Note  that  the  v vectors  corresponding  to  a characteristic 
value  X of  multiplicity  greater  than  1 may  be  assumed  to 
have  a special  form,  at  least  if  the  coordinate  vectors  are  appro- 
priately ordered.  For  definiteness,  suppose  X^  is  of  multiplicity  2. 

Let  v^  and  v^  be  linearly  independent  modal  vectors  corresponding  to  A ^ . 
Then  any  linear  combination  of  v^  and  is  also  a modal  vector 


r 


corresponding  to  A^.  From  this  fact  it  follows  that  there  are 
linearly  independent  vectors  and  v2  corresponding  to  A^  such 
that,  see  Fq.  (A15),  the  first  component  of  is  not  zero,  the 
second  component  is  zero,  the  first  component  of  is  zero, and 
the  second  component  is  not  zero.  Corresponding  to  these  vectors  v^ 
and  of  special  form  there  are  vectors  Uj  and  which  satisfy  the 
orthogonality  and  normality  conditions,  Eqs.  (A17)  and  (A18). 

Thus  whatever  the  degree  of  multiplicity  of  A greater  than  1, 
we  may  assume  that  there  is  a linearly  independent  set  of  v vectors, 
the  v portion  of  which  has  the  special  form  just  described  for  some 
consecutive  set  of  components.  Moreover,  corresponding  to  these  v 
vectors  of  special  form  there  are  u vectors  which  satisfy  the  orthogonality 
conditions  Eq.  (A17).  The  normality  conditions  can  be  met  by  appro- 
priately adjusting  (there  are  infinitely  many  ways)  the  magnitudes  of  the 
u and  v vectors. 

Let  A * A (t)  denote  the  diagonal  matrix  of  order  n for  which  the 


ktn  diagonal  element  is  exp(A^t).  Set  (see  Fq.  (A19)) 

U(t)  = U A 


U(t)  * U AA  = U AA 


(A21) 


(A22) 


AU  ♦ BU  • [AU  A ♦ BU]  A - 0 (A23) 

We  now  seek  a particular  solution  y(t)  of  Fq.  (A9)  where  we 


suppose 


y(t)  - u(t)z(t) 


(A24) 


and  the  vector  function  z(t)  is  to  be  determined.  Substitute  y(t)  and 
y(t),  as  determined  from  Fq.  (A24) , into  Fq.  (A9)  and  use  the  results, 
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Eq.  (A23)  to  obtain 

AU  A z(t)  = £(t) 

Multiplying  Eq.  (A25)  by  A**(t)V^  we  obtain 


t ~ 
A" 

0 


Hence 


and 


z(t)  = | 

y(t)  = 

Partition  the 

matrices 

U = 

U11  U12_ 
U21  U22 

VT  = 

^11  V12 
V21  V22 

U A(t  - x)\rf  (x)  dT 

0 


(A25) 


(A26) 


(A27) 


, A(t  - t)  = 


rT 


All  0 

0 A 


22 


(A28) 


This  partitioning  is  in  conformity  with  the  partitioning  of  the 
vectors  u,  v and  y,  Eqs.  (All),  (A15)  and  (A8).  One  readily  obtains 
from  Eq.  (A27) 

XjM  - f [UjjAjjCt  - t)v12  * u,^2(t  - T)V22]fW  dT 

(A29) 

It  is  clear  that  the  term  U-^A^ft  - T)vi2^T)  represents 


[exp  Aj(t  - t)u1  •••exp  Am(t  - x)um] 


Vj  f(x) 


(A30) 


= [exp  Xx Ct  - x)UlvJ  •••  exp  Ajt  - x^vj]  f(x) 


A similar  expression  is  obtained  for  the  term  22 (t  ' T)V22f(T). 
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Thus  we  obtain,  in  dyadic  form,  the  expression  for  x^(t) 


n rt 
= E 

k=l  J0 


U]Vk  exp  Ak(t  - x)f(x)  dr 


(A31) 


Note  that  the  subscript  on  x has  changed  from  1 to  0.  The 

use  of  the  subscript  1 is  discontinued  because  we  no  longer  care 

to  emphasize  the  relation  to  Eq.  (A8).  The  subscript  0 is  used  because 

xQ(t)  as  given  by  Eq.  (A31)  is  the  particular  solution  which 

satisfies  Eq.  (Al)  with  zero  initial  conditions. 

Take  f(t)  = h(t)r^  where  h(t)  = 0 for  t < 0,  h(t)  = 1 for 

t _>  0 and  r^  denotes  the  jth  coordinate  vector.  The  jth  component 

of  r.  is  1 and  all  other  components  are  zero.  Also,  let  v-t  denote 
J JK 

the  jth  component  of  the  vector  v^.  From  Eq.  (A31)  the  step  response 

HE  (t)  due  to  a unit  step  excitation  at  the  jth  station  is 
t 


Hj  (t)  = | ^ U]Vjk  exp  Ak(t  - t)  dx 


(A32) 


That  is 


n 


Hj(t) 


e,  ci/y  vjk  * £ Vjk exp  (V° 


k*l 


(A33) 


The  impulse  response  Ij  (t)  is  obtained  from  the  step  response  FT  (t)  by 
differentiation.  Thus  differentiating  Eq.  (A33)  gives 

:j  W * Vjk  exp  (V1 


(A34) 
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The  response  x^(t)  to  the  harmonic  excitation  r exp  (iait)  is 


XW(t)  = ■Jj.j  lVk  • r CXP(Xkt}| 


(A35) 


* ,"kvk  * r cxp(iwt)]/(ii»  - Ak) 


In  Eq.  (A35)  the  response  x (t)  is  expressed  as  the  "sum"  of  two  sums. 
The  first  sum  is  a linear  combination  of  solutions  of  t he  homogeneous 
equation,  Eq.  (A2),  called  the  complementary  function.  The  second  sum 
is  a particular  solution  of  Eq.  (A l ) with  f(t)  = r exp  (iwt).  If  the 
characteristic  values  A^  are  complex  numbers  with  Re  | A ^ | < 0 then  the 
complementary  function  goes  to  zero  as  t becomes  large.  ’Hie  complemen- 
tary function  is  a transient  and  the  particular  solution  is  called  the 
steady  state  solution. 

Whether  the  complementary  function  is  transient  or  not,  set 


M ... 

y = |j:  u. v,./ ( im  - A.)l  • r 


(A36) 


Then  y exp(iwt)  is  a particular  solution  of  Eq.  (Al)  when  f(t)  - r exp(iuit). 
For  the  time  being  we  will  refer  to  y exp(imt)  as  the  steady  state  solution. 
Many  procedures  for  determining  values  for  structural  parameters  are  based 


on  Eq.  (A 36).  It  follows  that 

(-<o^M  + K + iw  Cl  y = r 
and  from  this  equation  that 

[-rn^M  ♦ K - iw  C]  y ® r 
Then  (by  linearity) 


(A37) 


(A38) 


(Mir  ♦ CP  ♦ K|  (y  exp(iiot)  - y exp  (-imt)l/2i  = r sin  wt 

(A39) 
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That  is,  [y  exp(icot)  - y exp(-iwt)]/2i  is  a particular  response  to  the 
sinusoidal  excitation  r sinwt. 

TCiis  response  to  the  sinusoidal  excitation  can  be  rewritten  as 
[y  exp(iwt)  - y exp(-iu>t)]/2i  * [ (y  - y)/2i]  cos  wt  + [ (y  + y)/2]  sin  wt 

■ Im  [y]  cos  tut  + Re  [y]  sin  wt  (A40) 
From  this  equation  it  is  clear  that  if  the  steady  state  response  to  a 
sinusoidal  excitation  r sin  o>t  is  known  then  the  steady  state  response 
to  the  harmonic  excitation  r exp(iwt)  is  known  and  conversely. 

Let  denote  the  j—  component  of  Re  [y]  and  bj  the  j—  component 
of  Im  [y],  Then  in  the  usual  fashion  we  have 

a.  sin  wt  + bj  cos  wt  = [a2  + b2]1^2  sin  (iot  + 6j)  (A41) 


tan  0^  = bj/aj 


(M2) 


Eq.  (A41)  provides  an  alternate  expression  for  the  steady  state  response 
to  a sinusoidal  excitation. 

The  Eqs.  (A36)  and  (A37)  contain  an  important  relationship,  namely 


[-w2!*  + K + iw  C]  ^ = E ujcv£/ 


(A43) 


The  right  hand  side  of  Eq.  (A43)  is  the  frequency  response  function. 

If  the  modal  vectors  u^  and  v^  and  the  characteristic  values  X^  are 
known,  then  Eq.  (A43)  can  be  used  to  determine  the  matrices  M,  C and  K. 
We  must  remember  that  Eq.  (A43)  is  not  true  for  arbitrary  modal  vectors 
u^  and  v^.  The  Eq.  (A43)  was  obtained  by  assuming  the  vectors  u^  and  v^ 
appropriately  normalized. 

Let  f j (t)  denote  the  j—  component  of  the  vector  f(t)  and  as  before, 


40 


Vjk  the  jth  component  the  vector  v^.  Then 


v‘f(t)  = E v f (t) 
K j*i  3K  3 


(A44) 


Using  this  result  we  can  rewrite  Eq.  (A31)  as 
n m 

xQ (t)  * J E uR(  E Vjkfj(x))  exp  Xk(t.  - x)dx 
n ^ 3=1 


r m n 

J [E  (E  u^  exp  Ak(t  - x))f.(x)]dx 

n 3=1  ^=1 


That  is,  we  obtain 


f m 

xQ(t)  = j (T  Ij(t  - x)  • fj(x))dx 

0 j-1 


(A4S) 


The  Eq.  (A45)  expresses  the  response  xQ(t)  as  a sum  over  j of  the 
convolution  of  the  jth  impulse  response  with  the  jth  component  of 


the  excitation. 
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APPENDIX  B 

THE  FOURIER  INTEGRAL,  FOURIER  SERIES  AND  TRIGONOMETRIC 
INTERPOLATING  POLYNOMIALS 

In  this  section  we  exhibit  the  relation  between  the  Fourier 
series  representation,  a trigonometric  interpolating  polynomial 
representation  and  the  Fourier  Transform  of  a real  valued  function 
f(t).  The  results  of  this  section  may  be  found  or  ferreted  out  of 
standard  references  such  as  [7,  15,  16].  They  are  included  here  for 
completeness  and  ready  reference.  Let  us  assume  that  the  function 
f(t)  satisfies  whatever  conditions  needed  so  that  any  indicated  operations 
on  or  representations  of  f(t)  are  valid. 

If  f(t)  is  periodic  of  period  t then  f(t)  is  representable  as 
a trigonometric  series. 

00 

f(t)  = aQ/2  ♦ Z (an  cos  nut  + bn  sin  nut)  (Bl) 

where  u = 2tt/t.  The  Fourier  coefficients  are  given  by  the  equations 

a = (2/t)  ( f(t)  cos  n u dt 

n Jo 

b = (2/t)  [ f(t)  sin  n u dt  (B2) 

n JO 

Substituting  the  complex  form 

cos  nut  ■ [exp  (inut)  ♦ exp  (-inut)]/2  (B3) 

and 

sin  nut  * [exp  (inut)  - exp  (-inut)]/2i  (B4) 


j 
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in  Hq.  (Bl)  and  collecting  like  terms  in  exp  (inu>t)  and  exp  (-inwt) 
one  obtains 

00 

fit)  - aQ/2  ♦ (1/2)  E^  [(an-ibn)  exp  (inwt)  ♦ (an*ibn)  exp  (-inwt)] 


Next,  setting 


C0  ■ a</2 

cn  ' <an  - ib„>'2 

C-„  - (an  * ib„>'2 


(note  C = C ) one  obtains 
-n  n 


f(t)  = E Cj.  exp  (ikcot) 


the  Fourier  series  representation  of  f(t)  in  complex  form.  The  complex 
Fourier  coefficients  are  given  by 

C.  ■ — [ f(t)  exp  (-inait)  dt  (B8) 

K T J 0 

A 

A function  f(t)  defined  by  the  equation 

« N „ N * 

f(t)  C,  exp  (iukt)-kJ1  exp  (2iukt/x)  (B9) 

is  a trigonometric  interpolation  polynomial  in  conplex  form.  The  complex 
coefficients  Ck  are  determined  by  the  fimction  values  assigned  to  f(t)  at 
the  points  t a jx/N  for  j * 1 to  N.  First,  we  obtain  a rule  for 

A 

confuting  the  coefficients  C^.  Next,  if  f(t)  is  a periodic  function 

A 

of  period  t and  if  f(t)  is  the  trigonometric  interpolating  polynomial 


satisfying  the  conditions 


f ( j t/N)  - f (j t/N) 


(BIO) 


! 


for  j * 1 to  N,  we  obtain  the  relation  between  the  coefficients 
of  the  interpolating  polynomial  and  of  the  Fourier  series. 

The  desired  results  are  obtained  from  readily  established  properties 
of  the  N,  roots  of  unity.  The  principal  N-—  root  of  unity  is 


r^  = exp  (2iri/N) 


(B11) 


(B12) 


r J - 1 = 0 


+ rx  + 1) 


(ri  - DCrJ-1  ♦ r?'2  ♦ 


Since  for  N / 1,  r^  j*  1 it  follows  that 


rN-l  + rN-2  + ...  + r + i . 0 


rj,  = r^  = exp  (2irik/N). 


for  k = 1 to  N (or  k = 0 to  N - 1,  as  convenient).  Then 


and  it  follows 


r?1'1  + i^'2  + •••  + r,  + 1 = 0 
k k k 


(B13) 


(B14) 


(B15) 


(B16) 


IU*.« 


1^35 


'"-“‘TaF 


provided  k N (or  0)  and  for  k *=  N (or  0) 
N-l  N-2 

rN  * rN  + - ♦ rN  * 1 a N 


(B17) 


From  F.qs.  (B15)  and  (B14)  it  follows  that 


r > r,  + ...  ♦ rM  = 0 


(B18) 


and  from  F.q.  (15) 


i k 

rj  = r- 
k 3 


(B19) 


Note  also  that 


Tj.  = exp  (-2irik/N)  = exp  (2niN/N).exp  (-2trik/N) 
= rN-k 


(B20) 


Consider  the  matrix 


12  3 

r r v 
i i i 


M 


12  3 N 

r2  r2  r2  — r2 


r1  r2  r3 
rN  N N 


N 


(B21) 


Since  r£  = r^  the  matrix  M is  symmetric.  The  matrix  product 


M-M  = [sj£] 


(B22) 


where 


N -k  t 

sie  a 1 ri  * Tl 

k=1  J k 


(B23) 
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1)  If  4 - j then 
N 


s..  * E exp  (-2irijk/N)  exp  (2irijk/N)  = N 
k*l 


(B24) 


2)  If  l + j then 
N 


rN-j 


k N k 
Tl  rN+t-j  = 0 


k-1  " J * lc-1 
since  rN+£_j  f 1.  It  follows  from  Eqns.  (B24)  and  (B25)  that 

MM-MM-NI 

From  Eq.  (B9)  for  t = jr/N  we  have 

A N A 1, 


(B25) 


(B26) 


f(jx/N)  - £ C.  r* 
k=l  K J 


(B27) 


for  j = 1 to  N.  This  system  of  equations  can  be  expressed  in  vector  form 

f (t/N) 
f (2t/N) 

I 

= M 


f(T) 


It  follows  from  Eq.  (B26)  that 


r*  1 

C 


, a 


Sj 


(1/N)M 


C “I 
f (t/N) 


f (t) 


(B29) 
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Eqns.  (B28)  and  (B29)  exhibit  the  relation  between  the  coefficients 

A A 

and  the  values  of  the  function  f(t)  at  the  points  t = jx/N  for  J*  1 
to  N.  Thus  if  we  have  any  (real  valued)  function  f(t)  which  is  well  de 

A 

fined  at  the  points  jx/N  for  j = 1 to  N and  if  we  set  f(jx/N)  = f(jx/N) 

A 

then  Eq.  (B29)  determines  the  coefficients  and  Eq.  (B9)  defines  a 

A 

periodic  function  of  period  x which  satisfies  the  conditions  f(jx/N)  = 
f(jx/N)  for  j = 1 to  N. 

From  Eq.  (B29)  we  see  that  the  coefficients  are  obtained  as 

the  product  of  a vector  multiplied  by  a matrix.  If  one  defines  an 

operation  as  a multiplication  of  two  (complex)  numbers  followed  by 

an  addition,  then  the  product  of  an  N-vector  by  a matrix  of  order  N, 

2 

in  general,  requires  N operations.  An  algorithm  is  developed  in  [17] 
which,  in  principle,  determines  the  C^'s  in  less  than  2 N log2  N operations. 
Inplementations  of  this  algorithm  are  referred  to  as  Fast  Fourier 
Transformers. 

Let  f(t)  denote  a periodic  function  of  period  x and  let  f(t) 
be  the  trigonometric  polynomial  which  interpolates  f(t)  at  the  points 
t = jx/N  for  j = 1 to  N.  For  t = jx/N,  for  j = 1 to  N we  have  by  Eq.  (B7) 


' f(jx/N)  = £ exp  (27rikj/N) 

-00 

00 

= E C,  [exp  2itik/N]:) 

-00 

Now  we  can  express  the  integer  k,  -°°  < k < as 
k = l • N + k. 


where  0 < kj  < N - 1.  Hence 


(B30) 


exp  27rik/N  * exp  2irikj/N 


(B31) 


(B32) 
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It  follows  that 


E CjJexp  2irik/N] 


Let  f(t)  be  defined  only  for  0 < t < t.  Here  now  we  will 
exhibit  the  relation  between  the  Fourier  coefficients  of  some  periodic 
extensions  of  f(t)  and  particular  values  of  the  Fourier  transform 
of  f(t).  First,  however,  let  us  describe  what  we  mean  when  f(t)  is 
extended  periodically  "in  the  standard  way”.  It  is  clear  that  for 
any  value  of  t,  -<*>  < t < 00  we  can  write 


The  complex  Fourier  coefficients  of  f(t)  are  given  by  Eq.  (B8) 


(It  is  instructive  to  plot,  or  at  least,  draw  some  vertical  line 


1 1 ■■ 


segments  representing  t|C^|  at  points  kw  on  a horizontal  axis  for 
k = 0,  ♦ 1,  ♦ 2,  •••. 

Suppose  we  define 


(B38) 


fj(t)  * f(t)  for  0 < t < t 
fj(t)  = 0 for  t < t < 2t 

and  then  extend  f^(t)  periodically  in  the  standard  way.  The  complex 
Fourier  Coefficients  for  f^(t)  are 

2t  = | f(t)  exp  (-ikiot/2)  dt 


(B39) 


Observe  that 


c»  - ** 


(B40) 


for  k = 0,  +1,  +2,  •••.  We  note  that  doubling  the  interval  from 

t to  2t  has  halved  the  distance  between  successive  plots  of  the  values 

.(1) 


2t|C£  J|  as  compared  to  successive  plots  of  the  values  t | Cfc | . Repeating 
this  process,  that  is,  increase  the  interval  from  2t  to  4t  define 


f2(t)  = f(t)  for  0 < t < t 
* 0 for  t < t < 4t 


(B41) 


and  extend  f2(t)  periodically  in  the  standard  way  then  the  distance 
between  successive  plots  of  the  values  4t|c£2^|  is  half  the  distance 
between  successive  plots  of  the  values  2t  |C^  | . Thas  we  can  make 
the  number  of  complex  Fourier  Coefficients  which  plot  in  some  fixed 
interval  about  the  origin  as  large  as  we  please  by  choosing  the  interval 
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r 


r 


-r^zssrz 


t < t < b sufficiently  large. 

If  now  we  define  the  Fourier  Transform  of  f(t)  by  the  equation 


- f f(t)  e’iwt  dt 
J -OO 


f f(t)  e'iu)t  dt 
In 


then  it  is  clear  that 
4>  C2tt1c/ t ) - tC^ 

for  n ■ 0,  ♦ 1,  ♦ 2,  •••  and  that 


(B42) 


(B4S) 


♦ (ifk/x)  - 2tC£ 


(1) 


(B44) 


for  n ■ 0,  + 1,  ♦ 2,  •••,  and  so  on. 

The  Fourier  series  representation  and  Fourier  transform  of  a function 
f(t)  are  not  readily  available  for  use  in  general  because  of  the 
integration  process  required  for  the  determination  of  the  Fourier 
coefficient  and  the  Fourier  transform.  However,  in  view  of  the 
relations  exhibited  above  between  the  Fourier  coefficients,  certain 
values  of  the  Fourier  transform  and  the  coefficients  of  the  trigonometric 
interpolating  polynomial  it  is  clear  that  if  the  trigonometric 
polynomial  interpolates  f(t)  (or  an  appropriate  periodic  extension 
of  f(t))  at  sufficiently  many  points,  the  coefficients  of  the  trigonometric 
polynomial  will  be  good  approximations  of  the  Fourier  coefficients  and, 
when  properly  scaled,  of  these  certain  values  of  the  Fourier  transform 


; 1 


of  f(t). 
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APPENDIX  C 

SYSTEMS  OF  LINEAR  DIFFERENCE  EQUATIONS 
In  this  section  we  review  the  material  which  is  basic  for  the 
two  methods  for  determining  certain  structural  parameters  given  in 
[8,9]  and  [10]  respectively. 

For  the  first  method  given  in  [8,9]  the  system  of  m second  order 
equations 

Mx  + Cx  = Kx  = f (Cl) 

is  exchanged  for  the  system  of  n = 2m  first  order  equations 


x2  * -M'1^  - M_1Cx2  + M_1f 
This  system  may  be  written  as 
y = Ay  + f 


(C2) 


(C3) 


where  the  matrix  A in  partitioned  form  is  given  by 


A 


0 I 
-M_1K  -M^C 


(C4) 


Suppose  we  know  a nontrivial  solution  y(t)  of  the  homogeneous 
equation 


y - Ay 


(C5) 


For  n different  values  of  time  t^,  k • 1 to  n set 

yk  * y(V 


(C6) 
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and 


yk  - dy/dt|t  = (C7) 

From  Eq.  (C5)  we  obtain  the  matrix  equation 

txi  ynl  ■ Hy1  •••  ynl  (C8) 

It  follows  that  if  the  matrix  [y^  •••  ynJ  is  nonsingular  then 

A - [yx  •••  yn]  [yx  •••  y^'1  (C9) 

We  shall  see  that  the  homogeneous  solution  y(t)  must  be  a linear 
combination  of  the  n linearly  independent  solutions  defined  by  the 
modal  vectors  and  characteristic  values.  Appendix  A.  Thus  we  may  write 

y(t)  = ux  exp  (Axt)  + •••  + un  exp  (Xnt)  (CIO) 

where  here  denotes  either  the  k—  modal  vector  or  the  zero  vector. 
For  brevity  let  us  write  also 

zjk  = exp  (XjV  (Cll) 


Then  the  matrix  [y^  • • > 
matrices 

fyj  — y„J 


y ] may  be  expressed  as  the  product  of  two 


fQl  fl2 


u ] 
nJ 


‘11  12 


• • 4 


nl  n2 


In 


“nn 


(C12) 
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It  is  clear  from  Eq.  (C12)  that  if  for  any  values  of  k,  is  the 
zero  vector  then  the  matrix  [y^  •••  y ] is  singular. 

The  next  task  is  the  determination  of  functions  y(t)  and  y(t) 
which  satisfy  the  condition  y = Ay  when  the  given  information  is  the 
second  derivative  x(t)  of  a function  x(t)  which  satisfies  Mx  + Cx  ♦ Kx  = 0. 
Let  us  set 

x(t)  = $(t)  + (Cl 3) 

Hence  by  Eq.  (C2)  we  have 


and 


x2(t)  = S(t)  + *0 

ft 


XjCt)  = x2(t)  = 


j [SCO  + SQ]  dx  + Cj 
0 


(C14) 
(Cl  5) 


Xj (t)  = | dxj  [X(s)  + ds  + Cjt  + C2 
0 0 

The  quantities  XjCt),  Xj(t),  x2(t)  and  x2(t)  defined  by  Eqs. 
(C15)  and  (C16)  satisfy  the  equation 


(C16) 
(C14) , 


Xjft) 

*2Ct). 


- A 


Xj(t) 

x2Ct) 


(C17) 


Now  since  Eq.  (C17)  holds  in  particular  for  t * 0,  the  Eq.  (C17) 


simplifies  to 


■ 

t 

L fT  1 

l*(0  ♦ S0]  dT 

Jo 

■A 

lodTJ0^(s)  + *0]  ^ + Clt 

*(t) 

f (S(T)  . *„]  dT 

- 

u 

(C18) 


All  the  quantities  in  Eq.  (C18)  are  known  except  the  integration 
constant  Cj.  Replace  t by  at  in  Eq.  (C18).  One  then  has 
rat 


|o  [*(t)  * *Q]  dr 
X(at) 


r*r 

Jo  Jo 


[£(s)  + Xq]  ds  + Cjat 


fat[X(T) 
J n 


+ dT 


CC19) 


Finally  multiply  Eq.  (C18)  by  a and  substract  from  Eq.  (C19) . 

One  obtains  an  equation  which  contains  the  matrix  A but  no  other 
unknown  quantities. 

Thus  from  the  above  discussion  we  obtain  for  the  vector  functions 


y(t)  and  y(t) 


y(t) 


| [*(t)  + SQ]  dr  - aj^[*(x) 

X(at)  - X(t) 


+ XQ]  dx 


(C20) 


and 


y(t) 


£ dxj*[$(s)  + ^01  ds  - aj*dxj*[*(s)  + SQ]  ds 
j01  [SCO  + XQ]  dx  - a | [X(x)  + XQ]  dx 


(C21) 


The  quantities  y(t)  and  y(t)  are  computed  from  the  measured  quantity 
x(t).  They  are  obtained  by  a numerical  integration  process  which 
has  been  made  somewhat  more  laborius  by  the  necessary  elimination  of 
an  integration  constant.  Convenient  values  for  a are  readily  apparent. 
Let  us  surmarize  the  discussion  of  this  section  up  to  tnis  point. 

We  noted  above  that  the  homogeneous  system  of  equations  (f  = 0),  Eq.  (Cl) 


has  a set  of  linearly  independent  solutions.  Next  given  a linear 
combination  of  all  these  linearly  independent  solutions,  the  matrix  A 
of  an  equivalent  first  order  system,  Eq.  (C5),  can  be  determined  from 
the  known  values  of  this  solution  of  the  homogeneous  system  Eq.  (Cl). 

The  second  method,  given  in  [10],  is  similar  in  principle  to  the 
first.  However,  it  avoids  the  numerical  integration,  Eqs.  (C20)  and 
(C21)  required  by  the  first  method.  The  procedure  which  enables  one 
to  bypass  the  nunerical  integration  is  described  in  [ 7,  pp.  272-280]. 

This  second  method  is  based  on  the  fact  that  corresponding  to  the 
homogeneous  differential  equation,  Eq.  (Cl),  there  is  a conpanion 
difference  equation  (of  second  differences)  having  the  same  solution 
set.  Then,  as  above,  given  a linear  combination  of  all  the  members  of  a 
linearly  independent  set  of  solutions  of  the  homogeneous  equation 
(f  = 0),  Eq.  (Cl),  one  can  determine  an  equivalent  equation  of  first 
differences.  Let  us  explain  in  detail  the  remark  just  made. 

Consider  a system  of  equations  of  second  differences 

x (t  + 2h)  ♦ Bx(t  ♦ h)  + Cx(t)  - 0 (C22) 

where  B and  C denote  matrices  of  order  m and  x(t)  is  a vector  valued 
function  of  t of  dimension  m.  It  is  clear  that 

x(t)  - u exp  (Xt)  (C23) 

satisfies  Eq.  (C22)  if 

n * exp  (Xh) 
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satisfies  the  characteristic  equation 


det[n  I ♦ Bn  + C]  - 0 


(C24) 


and  u is  a corresponding  characteristic  vector  satisfying  the  condition 


[n  I ♦ Bn  ♦ C]u  * 0 


(C25) 


Xj(t  ♦ h)  - x2(t)  ■ 0 

x2(t  + h)  + CXj(t)  + Bx2(t)  = 0 

The  system  of  equations,  Eq.  (C26)  can  be  written  as 


(C26) 


xx(t  + h) 
x2(t  + h) 


1 


B x2(t) 


(C27) 


a system  of  equations  of  first  differences  which  is  equivalent  to 
Eq.  (C22). 


x(t)  * I Ul  exp  (X.t)  ( 

k=l  K K 

n ■ 2m,  be  a sun  of  linearly  independent  functions  which  satisfy 
the  conditions  Eqs.  (C24)  and  (C25).  If  the  values  of  x(t)  are 
known  at  times  t^  and  t^  + h for  k * 1 to  N then  we  have  the  matrix 
equation 


(C28) 


Xl^l  + “*  Xl^n  + h)  0 

x2(tj  + h)  •••  x2(tn  + h)  C 


X1(V  •"  xl(tn> 

x2(tj)  •••  x2(tn) 


(C29) 
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: 1 - . 1 ■■■' 


Eq.  (C29)  is  solvable  for  the  matrix  of  the  system  (C27)  if  the 
matrix 


*l(‘)  •••  Vtn> 

x2(t)  — Xj(t^) 

► - 


is  nonsingular 

From  this  discussion  it  is  clear  that  given  a function  x(t)  as 
defined  by  Eq.  (C28)  one  can  determine  a system  of  equations  of  first 
differences  which  the  n- dimensional  vector  valued  function 


y(t) 


x(t) 

x(t  + h) 


satisfies.  Thus  we  determine  the  matrices 


(C30) 


$ = 


and 


x(t,)  •••  x(tn) 

x(tx  + h)  •••  x(tn  + h) 


(C31) 


$ 


x(tj  + h)  •••  x(tn+h) 

x(tl  + V ***  x^n  + 2h) 


(C32) 


Then,  if  4>  is  nonsingular,  y(t)  defined  by  Eq.  (C30)  satisfies  the 
system  of  equations  of  first  differences 


y(t  ♦ k)  - * 4>’1  y(t)  * 0 


(C33) 
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APPENDIX  D 
COMPUTER  PROGRAM 


We  now  describe  the  organization  of,  and  the  calling  procedure 
for  invoking  the  use  of,  a collection  of  subprograms,  written  in 
FORTRAN  to  implement  the  method  described  in  Section  II.  It  accepts, 
as  input,  the  frequencies  and  the  vectors  and 
(obtained  from  experimental  data),  together  with  various  other  data, 
and  produces,  as  output,  the  characteristic  values,  x^,  the  modal 
vectors,  and  V^,  and  the  mass,  damping,  and  stiffness  matrices, 

M,  C,  and  K.  A section  on  program  testing  is  also  included. 

ORGANIZATION  - Subroutine  AAS,  called  by  the  user,  allocates 
virtual  storage  in  the  work  array,  W,  and  calls  subprogram  BBS. 

BBS  monitors  and  controls  the  computations  and  calls  into  use  the 
following  subprograms: 

EIGEN  - computes  an  approximation  to  a characteristic  value, 

X^,  given  the  associated  frequencies  and  the  vectors 

and  Y^2  (or  and  Y^).  See  Eqs.  (7),  (8),  (14)  and  (15),  Section 
II. 

VUS  - computes  an  approximation  to  a modal  vector  and  the 
first  component  of  the  associated  vector  V^,  given  the  frequencies 

and  ^2 » the  characteristic  value,  x^,  and  the  vectors  Ykl  and 
Yk2  (or  Ykl  31141  Yk2)#  866  Eqs*  (10).  (ID, (16),  and  (17),  Section  II. 

YVECS  - computes  Ykl  and  Yj^,  given  the  vectors  Y^  and  Y^7, 
the  frequencies  u^,  and  the  current  approximations  to  the 
characteristic  value,  x^,  and  the  modal  vectors  and  V,..  See 
Eq.  (13),  Section  II. 


Iddl 
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RECVR  - computes  the  system  matrices  M,  C,  and  K,  given  the 
characteristic  values,  and  the  modal  vectors  and  Vj,.  See  Eqs. 
(38) -(43),  Section  II. 

INV,  MTMT,  LU3  - used  by  RECVR  to  perform  various  matrix 
operations. 

Together,  these  subprograms  utilize  10  cells  of  temporary  work 
space  in  blank  COM  ION,  in  addition  to  the  work  array,  W. 

GVLLING  PROCEDURE  - Die  calling  program,  written  by  the  user, 
will  contain  the  equivalent  of  the  following  FORTRAN  statements: 

EXTERNAL  DATA 

COMPLEX  LAMBDA  (M)  ,U(M,M)  ,V(M,M) 

REAL  EM  (M  ,M)  , C (M  ,M)  , KA  (M  ,M)  , W (NW) 


where : 


CALL  A/VS  (M,DATA,TOL,IMAX, LAMBDA, U,V,EM,C,KA, I, W) 

M »s  an  integer  input  variable,  indicating  an  m by  m sys- 
tem of  second  order  differential  equations. 

DATA  is  the  name  of  a user  written  subroutine  subprogram,  called 
by  BBS  to  obtain  certain  input  data.  The  form  of  this  is 
described  in  a later  section.  The  name  "DATA"  is  symbolic, 
of  course,  and  can  be  any  legal  subroutine  name  that  does 
not  conflict  with  ones  in  use  in  the  collection.  (See 
list  under  ORGANIZATION.) 

TOL  is  a real  input  variable,  the  tolerance  used  to  control 
program  convergence.  Let  be  the  ith  iterate  of  x^. 

TTie  process  is  said  to  converge,  no  further  iterations 
being  performed,  if  1^  - 1 | < TOL,  for  all 

k,  k=l,2,...,m.  Some  care  should  be  exercised  here  to 
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ensure  that  TOL  is  not  unreasonably  small. 

IMAX  is  an  integer  input  variable,  the  maximum  nunber  of  itera- 
tions allowed.  15^IMAX<30  would  probably  be  reasonable. 

LAMBDA  is  a complex  output  array,  of  dimension  m,  containing  the 
characteristic  values,  k=l,2,..,m. 

U,V  are  complex,  m by  m,  output  arrays,  containing  the  modal 
vectors,  U^  and  V^,  respectively.  U(I,K)  will  contain  the 
I-th  component  of  vector  Uj.,  and  similarly  for  V(I,K)  and  V^. 

EM,C,KA  are  real,  m by  m,  output  arrays,  containing  the  mass  damp- 
ing, and  stiffness  matrices,  M,  C,  and  K,  respectively. 

I is  an  integer  output  variable,  containing  the  number  of 
iterations  required  for  convergence.  If  convergence  is 
not  achieved,  I is  set  to  IMAX  +1. 

W is  an  array  used  by  the  subprograms  as  temporary  work  space. 

2 

It  must  be  dimensioned  at  least  4m  +6m. 

DESCRIPTION  OF  SUBROUTINE  DATA  - this  subprogram,  written  by  the  user, 
is  called  by  BBS  to  provide  it  with  the  following  data: 

1.  the  frequencies  <aKl  and  k=l,2,...,m. 

2 . the  vectors  Ykl  and  Yj^.  lc=  1,2, .. . ,m;  and  of  length  m. 

These  are  obtained  from  experimental  data. 

This  subprogram  should  contain  the  equivalent  of  the  following  FORTRAN 
statements : 


SUBROUTINE  DATA  (M.CMEGA.Yl ,Y2) 
REAL  CMEGA(N) 

COMPLEX  Y1  (M,M)  ,Y2  (M,M) 
where: 


I 

I 

I 
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M is  an  integer  input  variable,  the  same  M as  provided  in  the 


initial  call  to  AAS. 


OMEGA  is  a real  output  variable,  of  dimension  N=2M,  which  will 
contain  and  u^-,,  as  follows: 

OMEGA(l)  will  contain 
0MEGA(2)  will  contain 


OMEGA (M)  will  contain  co^ 
OMEGA (M+l)  will  contain  u)j2 


OMEGA(N)  will  contain  to  - 

mZ 

Y1,Y2  are  complex,  m by  m,  output  arrays,  which  will  contain  the 
vectors  Ykl  and  Yk2,  respectively.  Y1(I,K)  will  contain  the 
I-th  component  of  vector  Ykp  I,k=l,2. . . ,m;  and  similarly 
for  Y2(I,K)  and  Yk2. 

PROGRAM  TESTING  - The  results  of  two  test  cases  with  m=4  are 
given  in  Tables  I and  n.  For  these  tests  the  values  of  Xk,  Uk  and 
Vk  were  specified.  Then  given  the  input  frequencies,  <0^  and  10^, 
we  computed  the  vectors  Ykl  and  Yk2,  using  Eq.  (2)  of  Section  II. 
Using  a tolerance  of  10  , convergence  was  achieved  in  11  iterations. 

The  following  tables  display  the  agreement  achieved  between  the 
specified  and  confuted  values  of  *k  and  Uj,.  Similar  agreement  was 
achieved  for  Vk,  obviously. 


' • 
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in  the  first  test  case  the  A,  are  all  different.  In  the 


second  case  Aj=A7.  The  remaining  A^  and  the  and  Vj.  are  the 
same  as  in  the  first  case.  The  specified  values  for  the  A,  and  U 


are  given  in  the  tables.  For  the  we  took  Vj  = i /T~ • 10 

V7  = i ST-  10_3U?,  \\  = i ST-  10_3U,  and  V,  = i/T-  10'\j„.  Note 


in  Table  2 that  the  computed  values  of  Uj  and  U7  are  linear  combina 
tions  of  the  specified  values  of  U.  and  U7. 


TABLE  I 


CHARACTERISTIC  VALUES  ALL  DIFFERENT 


Exact 


-0.2  + i 

-0.05  + 2i 


-0.1  + 3i 


-0.175  + 4i 


-13/9 


Computed 


-.20000000  + 1.0000000  i 
-.050000000  + 2.0000000  i 
-.10000000  + 3.0000000  i 
-.17500000  + 3.9999999  i 


.66666668  + .841  X 10'°  i 
1 .0000000  + .228  X 10'8  i 
-.173  X 10‘8  - .154  X 10'8  i 


-1.5000000  - .390  X 10  i 
.148  X 10'7  + .677  X 10'8  i 
1.0000001  + .177  X 10'8  i 


.66666665  - .313  X 10'8  i 
-1.4444444  - .214  X 10'8  i 
.137  X 10'7  + .123  X 10’7  i 


-.50000002  + .364  X 10'8  i 
.142  X 10'8  - .355  X 10'8  i 
-2.4999999  + .108  X 10'7  i 
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TABLE  II 

A CHARACTERISTIC  VALUE  OF  MJLTIPLICITY  2 


l 

I 

' 

\ 

I 

) 

» 

i 
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Data  Exact 


Computed 


U, 


U, 


U, 


-0.2  + i 
-0.2  + i 
-0.1  + 3i 
-0.175  ♦ 4i 
r 1 
2/3 

- 

1 

. 0 
1 

-3/2 

0 

1 

1 

2/3 
-13/9 
. 0 
’ 1 
-1/2 
0 

, -5/2 


-.20000000  + 1.0000000  i 
-.20000000  1.0000000  i 

-.10000000  + 3.0000000  i 
-.17500000  + 4.0000000  i 
1. 

-.52610556  - .704  X 10‘12  i 
.44948974  + .123  X 10_11  i 
.55051026  ♦ .285  X 10'12  i 
1. 

-2.7340926  ♦ .343  X 10'10  i 
-.56958119  - .253  X 10‘10  i 
1.5695812  - .232  X 10* 11  i 
1. 

.66666667  - .253  X 10* 9 i 
-1.4444444  + .401  X 10*9  i 
-.840  X 10"10  ♦ .117  X 10'9i 
1. 

-.50000000  - .374  X 10*11  i 
-.292  X 10*10  ♦ .631  X 10*10  i 
-2.5000000  + .427  X 10*9  i 


64 


1*13  MQE  IS  BfcST  QUALITY  PBACtlCAW* 

from  COPY  FURNISHED  TO  DDC '' 


SUBROUTI  NE  AAS  ( M,  D A T A,  TOL » I MX , L AMBOA  , U,  V,  FM,  C , KA  , I , W) 

C 

C THIS  SUBROUTINE  ALLOCATES  VIRTUAL  STORAGE  IN  THE  WORK 

C ARRAY,  W , DEPENDING  ON  THE  PAPAMETEP,  M,  AND  CALLS  THE  CORE 
C SUBROUTINE  BBS • 

C 

EXTERNAL  DATA 
01  HE  NS I ON  Hill 
K=2*m 
I Y 1=  1 ♦< 

IY2=I Y l^K 
IY=IY?*K 
K = K*M 
IYT=IY*< 

K=M*M 
IT3  = IY 
IT4=IT3*K 
IT5=IT4*K 

CALL  BBS  (H , P AT  A , T OL , I MAX, LAMBDA, U*V,EM,C»KA*I, W (1) 

1 ,W(IY1»,W(IY2),W(IY),W(IYT),H(IT3),M(IT4),W(ITE)) 

RETURN 

END 


f 


II 


<1 

ft 
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THIS  PAGE  is  BEST  QUALITY  PRACTICAL* 

Jgau  C0rY  EUKHISR^  r°  DD' 


SU3R0UT I NE  BBS  < M, DA T A , TOL , IMA  X , LAMBDA , U , V, EM, C , KA ,LL 
1 , OMEGA,  Y1,Y2,Y,YT,TT,TA,T5> 

THIS  SUBROUTINE  MONITORS  AMD  CONTROLS  ALL  C ALCUL  AT  I CAS , 
INC  LUOI MG  THE  READING  OF  THE  INPUT  DATA  VIA  THE  USER  WRITTEN 
SUBROUTINE  'DATA*. 

EXTERNAL  OATA 
DIMENSION  OMEGA(l) 

COMPLEX  T, LAMBDA (M) ,Yi(M),Y2<M),V<M,M),U<M,M)tY<M,M>,YT(M,M> 

READ  INPUT  DATA*  OMEGA'S  AND  RESPONSE  VECTORS  Y 

CALL  DATA  (M, OMEGA , Y , YT> 

FIRST  APPROXIMATION 

DO  20  K= 1 » M 

OMl=OMEGA(tO 

OM2=OMEGA(K+M> 

OBTAIN  EIGENVALUE,  LAMBDAOO 
CALL  EIGEN  (Y(i,K)  ,YT  ( 1,  K)  ,OMl,OM2,LAMBDA(K),M) 

OBTAIN  EIGENVECTORS,  U ANO  V 

CALL  VUS  <Y<i,K),YT(l,*0  , LAMBDAOO  ,OMi , OM2 , V < 1 , K ) ,U(1,K),M» 
20  CONTINUE 

REFINEMENT  PROCEDURE 

00  EE  L=1 »I MAX 
KN=0 

DO  AO  K=1,M 
OMl=OMEGA(K) 

0M2=0M£GACKfM) 

OBTAIN  NEW  Y VECTORS,  Y1  AND  Y2 
CALL  YVECS  <K , OMEQ  A , Y , YT  ,L  AMBDA  , V,U , Y1 , Y 2,  M ) 

OBTAIN  NEW  EIGENVALUE,  T = LAMBDAOO 


CALL  EIGEN  (Yl, Y2 , 0M1 ,OM2 ,T , M ) 


OBTAIN  NEW  EIGENVECTORS 


CALL  VUS  (Y1,Y2»T,0M1,0M2»V(1»K) , U ( 1 , K ) * Ml 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROM  COi  Y FURNlSHi&D  TO  UDC  


TEST  FOR  CONVERGENCE  AND  STORE  N^W  EIGENVALUE 
IF  ( CA DS ( LAMBDA ( K) -T I ,G£.  TOL)  KN=KN^1 

uc  lamboa<k>*t 

IF  CONVERGENCE  MAS  MET,  GO  TO  f<3 

IF  (KN  .EQ.  0>  GO  TO  o3 
SO  CONTINU" 

CONVERGENCE  FAILED 

LL  = INAXU 
GO  TO  70 

SET  ITERATION  COUNT  AND  COMPUTE  THE 
REMAINING  ELEMENTS  OF  VECTOR  V 

6C  LL=  L 

70  DO  83  K=1,M 

00  80  J - ? , M 

H V( J,K> =V(i,K)*U( J,K) 

RECOVER  MATRICES  M,  C,  < 

CALL  RECVR  01 , V ,LA  MOD  A,  M ,fM  , KA  , C , T 3 ,T  U ,T  5 > 

RETURN 

en  n 


SUBROUTINE  FI Gp  N < Y 1 , Y? , 0M1  , QM2  . T,  M) 


THIS  SUBROUTINE  COMPUTES  AN  APPROXIMATION  TO  LAMBOA(K), 
GIVEN  THE  FREQUENCIES  OMEGA ( K 1 ) ANO  0MEGA(K2»,  ANO  THE 
CURRENT  APPROXIMATION  TO  THE  VECTORS  Y(K1>  AND  Y <K?> . 

COMPLEX  Yl(  1»  , Y2(i)  , T,H  ,T2 
COMMON  T1,T2 
Ti  = (C.,0.) 

T2=<0.,0.» 

no  ic  1=1, m 

T1=T1*Y1III **2 
1C  T2=T2*Y2<I> **2 
T=  CSDPT (T1/T2I 
IF  ( A IMAG (T > ,LT,  0.  > T=-T 
T1=CM»LX(3. ,0Mi) 

T2 =CMPLX(0. ,0M2) 

T=  <T*T1-T2)/(T-1.I 

RETUFN 

ENO 
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SHIS  PAGE  iS  BEST  QUALITY  PRACIICA%J 

{BOM  Cui.  1 !■  Otu<i ISHEi)  TO  DDC  ^ 

SUBROUTINE  VUS  ( Y, YT , LAMBOA , OH1 ,OM2 ♦ V , U, M> 

THIS  SUBROUTINE  COMPUTES  VECTOR  U(K»  AND  THE  i-ST 
ELEMENT  OF  VECTOR  VltO,  GIVEN  THE  CURRENT  APPROXIMATION 
TO  THE  CHARACTERISTIC  VALUE  LAMBOA ( K > * THE  FREQUENCIES 
OMEGA ( <1 I AND  OMEGA <K2>,  ANO  THE  CURRENT  APPROXI MATION 
TO  THE  VECTORS  Y(K1>  ANO  Y(K2) . 

COMPLEX  Y(1>,YT<1> ,L AMQOA , V ,U< 1 > , 01 

COMMON  01 

01=YT<1I-Y(1» 

V=01/ C1./(CMPLX(0. ,0M 2) -LAMBDA) -1./ (CMPLX(0 . ,0M1) -LAMBOA) ) 

U( 1>=<1.,0. ) 

DO  10  L=  2»M 

UCL)=(VTtL> -Y<L> )/01 

CONTINUE 

RETURN 

ENC 


SUBROUTINE  YVECS  ( K , OMEG A , Y , YT , LAMBDA , V , U, Y 1 ♦ Y2 , M I 

THIS  SUBFOUTINE  COMPUTES  AN  APPROXIMATION  TO  THE 
VFCTORS  Y(K1>  ANO  Y(K2>,  GIVEN  THE  FREQUENCIES  OMEGA ( Kl) 

ANO  OMEGA < K2 ) » THE  RESPONSE  VECTORS  Y,  AND  THE  CURRENT 
VALUES  OF  ALL  LAMBOA  * S » V'S,  AND  U«S. 

DIMENSION  OMEGA(l) 

COMPLEX  Y(M,M> ,YT(M,M)  ,LAMBDA(1) ,V<M,M),U<M,M),Y1(1>,Y2<1>, 
1 T1,T2,S1,S2,T 

COMMON  T1,T2»S1,SZ»T 
T1=CMPLX(0. ,OMEGA<  K) ) 

T2=CMPLX(0.,OMEGA(K+M)) 

00  30  1=1, M 

SI =Y ( I , <) 

S2=YTCI,KI 
no  20  L = 1,M 
T=V(1,L)*U(I,L > 

IF  <L  . EO.  K)  GO  TO  10 
S1=S1-T/  (T1 -( LAMBDA  (L)  > ) 

S2  = S2-T/ (T2-< LAMBDA  <L>  > > 

S1  = S1-<C0NJG(T) )/< Tl-(CONJG( LAMBOA (L)  I )> 

S2=S2- ( CON  1G(TM/(T2-(C0NJG(  LAMBOA  (L))  >> 

CONTINUE 

Yl( I)=S1 

Y2(I»sS2 

CONTINUE 

PETURN 

ENO 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 

mM  COPY  FURNISHED  TO  UDC  - 

SUBPOUTINE  RECVR  ( U, V , LAMBDA, M , EM, «A, C ♦ Tl , T2, T3 ) 

THIS  SUBROUTINE  COMPUTES  THE  MASS,  DAMPING,  AND 
STIFFNESS  MATRICES,  M,  C,  AND  K,  GIVEN  THE  FINAL  VALUES 


C OF  ALL  LAMBDA'S,  V*S,  AND  U'S. 

C 

COMMON  T,L,VV,M2,N 

COMPLEX  UCM,M) ,V<M,M> ,LAMBOA(l» ,T,L,VV 
DIMENSION  FM(1) , KA ( 1 ) , C ( 1 ) , T 1 ( 1 ) , T 2( 1)  * T 3 ( 1 ) 
M2=M*M 

DO  10  1*1, M2 

10  Tl <I>=T2(I)=T3<I>=0. 

C 

C COMPUTE  MATRICES!  H<  0)  , -I*H»<0>, 

C 

DO  20  K = 1,M 
L= LAMBDA (K) 

N=0 

DO  20  J=1,M 
VV=V(J, K) 

DO  20  I = 1 ♦ M 
N=N  +1 

T=-U(I,K>*VV/L 
T1(N)=T1(N) ♦REAL(T) 

T = T/L 

T2 (N) =T  2 ( N)  ♦FEAUT) 

T=T/L 

T3 ( N) =T3(N> -REAL(T) 

20  CONTINUE 

OO  30  1=1, M2 
T1II)=2.*T1(I) 

T2J I) =2  , *T2 ( I ) 

T3 ( I) =2  »*T3  C I ) 

30  CONTINUE 
C 

C COMPUTE  REAL  MATRICES  M,  C,  K 

C 

CALL  INV  <T1,KA,M,EM) 

CALL  MTMT  ( M,  M ,M  ,T  2 , K A,  T 1) 

CALL  MTMT  ( M»M,M,K A, Tl,C I 
CALL  MTMT  < M,M ,M ,C ,T 2 , Tl ) 

00  40  1=1, M2 

40  cm*-cm 


CALL  MTMT  (M,M,M,KA,TT,T2) 
DO  SO  1=1, M2 
50  T3(I>*Tim  ♦T2(I) 

CALL  MTMT  C M,M ,M«T  3 , K A,EM) 

RETURN 

ENC 
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